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The  program  EROS  (Energy  Recovery  Optimization  System)  is  a  flow- 
sheeting  package  for  evaluating  and  optimizing  the  performance  of  simple 
networks  of  heat-exchangers.   Using  EROS  one  can  set  up  an  arbitrary 
structure  of  heat-exchangers,  stream  splitters  and  mixers.   Stream 
flow-rates  and  entry  and  exit  temperatures  may  be  specified,  free,  or 
bounded  above  and/or  below.   Phase  changes  may  be  allowed  to  occur. 
Requiring  no  more  user  input  (such  as  initial  guesses),  the  program 
gathers  together  the  modeling  equations  and  appropriate  inequality 
constraints.   It  then  develops  solution  procedures  repeatedly  in  the 
course  of  optimizing,  initially  to  aid  in  locating  a  feasible  point 
(which  need  not  be  specified  by  the  user)  and  in  the  final  stages  to 
take  advantage  of  tight  inequality  constraints  to  reduce  the  degree 
of  freedom.   The  program  is  written  in  standard  Fortran  and  is  cur- 
rently operable  on  IBM  360/67. 

The  development  of  this  program  is  relevant  from  the  point  of  view 
of  evaluating  and  optimizing  energy  recovery  systems  because  the  subject 
of  choosing  an  optimal  structure  has  recently  come  under  increasing 
attention.   A  tool  such  as  EROS  will  prove  very  useful  to  make  a  more 
thorough  analysis  of  the  candidate  structures  chosen  by  the  process  of 
synthesis  at  the  penultimate  and  final  stages.   Existing  energy  recovery 
schemes  may  be  reevaluated  with  the  intention  of  making  changes  to  make 
them  more  efficient.   EROS  can  also  perform  quick  reliability  studies 
accounting  for  fluctuations  in  stream  flow-rates  and  temperatures, 
a  common  problem  during  start-up  and  periodic  failures. 

As  is  the  case  in  the  optimization  of  most  chemical  engineering 
systems,  a  global  optimum  cannot  always  be  guaranteed.   An  attempt 


to  recognize  whether  the  discovered  optimum  is  a  global  one,  in  a 
large  system,  is  not  a  trivial  task.   Usually  the  problem  may  be 
side-stepped  by  making  several  optimization  runs  from  different 
starting  points  and  considering  the  best  result  as  the  required 
optimum.   A  systematic  procedure  is  presented  in  this  study  to 
investigate  global  optimality,  and  the  results  on  application  to 
simple  examples  prove  it  to  be  quick  and  effective. 

The  use  of  EROS  can  be  extended  to  perform  process  synthesis 
via  structural  parameters.   However,  the  presence  of  inequality 
constraints  inherent  in  the  system  can  easily  force  the  optimization 
to  stop  short  of  the  desired  result.   This  observation  has  not  been 
reported  before  and  attention  is  brought  to  it  in  this  dissertation. 


CHAPTER  I 
INTRODUCTION 

Most  of  the  existing  flowsheet ing  packages  for  chemical  engi- 
neering processes  are  based  on  Sequential  or  Simultaneous  Modular 
Approaches.   In  these  approaches  each  unit  is  modeled  by  writing  a 
computer  subroutine  which  converts  the  input  stream  and  equipment 
parameter  values  into  output  stream  values.   Systems  based  on 
Sequential  or  Simultaneous  Modular  Approaches  are  relatively  easy 
to  build  but  the  penalties  paid  are  the  lack  of  flexibility  in  the 
definition  of  the  problem  and  the  requirement  for  well-defined  user- 
specifications  . 

Equation  solving  approaches,  as  followed  in  EROS  (and  in  Leigh, 
Jackson  and  Sargent  (1974),  Hutchison  and  Shewchuk  (1974)  and  Kubicek, 
Hlavacek  and  Prochaska  (1976))  present  an  alternative  way  to  treat 
flowsheets.   The  flowsheet  is  represented  as  a  collection  of  non-linear 
equations  which  must  be  solved  simultaneously.   In  following  an  Equation 
Solving  Approach,  the  user  can  specify  many  of  the  values  for  both  unit 
inputs  and  outputs.   The  unit  equipment  parameters  are  then  calculated 
to  give  these  desired  transformations  of  inputs  to  outputs  by  the  unit; 
in  other  words,  the  unit  is  designed  to  meet  these  requirements.   Since 
EROS  also  contains  an  optimization  capability  the  rather  striking 
advantage  is  that  variables  for  whose  values  the  user  has  no  preference 
for  are  treated  as  the  degree  of  freedom  by  the  system.   EROS  incor- 
porates the  general  optimization  strategy  as  outlined  in  Westerberg 


and  deBrosse  (1973)  and  demonstrates  the  applicability  and  effec- 
tiveness of  their  algorithm. 

The  two  important  problems  in  the  design  of  energy  recovery 
systems  are  to  choose  the  configuration,  and,  given  a  configuration, 
to  choose  the  design  parameters  and  operating  variables.   A  recent 
review  on  the  effort  directed  to  choosing  a  configuration  can  be  found 
in  Nishida,  Liu  and  Lapidus  (1977).   In  choosing  a  suitable  configura- 
tion the  general  trend  has  been  to  evaluate  networks  using  the  heuristic 
of  setting  the  minimum  allowable  approach  temperature  to  20°F,  whereas 
the  economics  often  advocate  a  smaller  value.   In  addition,  when  a 
stream  is  split,  the  need  for  finding  the  optimal  value  for  the  split 
fraction  has  been  ignored.   Grossman   and  Sargent  (1977)  optimized 
several  heat-exchanger  networks  and  found  considerable  savings 
(sometimes  as  much  as  25%)  . 

The  problem  of  optimizing  a  heat-exchanger  network  to  obtain 
the  most  suitable  values  for  operating  variables  has  been  considered 
in  the  works  of  Westbrook  (1961),  Boas  (1963),  Fan  and  Wang  (1964), 
Bragin  (1966)  and  Avriel  and  Williams  (1971).   Typically  each  design 
problem  is  formulated  and  solved  for  as  an  optimization  problem. 
Many  investigators  (Hwa  (1965),  Takamatsu,  Hashimoto  and  Ohno  (1970), 
Henley  and  Williams  (1973),  and  Takamatsu  et  al.  (1976))  have  combined 
both  the  problems,  choosing  a  configuration  as  well  as  the  operating 
variables,  and  formulated  it  as  an  optimization  problem. 

All  the  methods  mentioned  for  optimizing  over  the  operating 
variables  require  the  problem  to  be  cast  into  a  mathematical  format. 
EROS  precludes  this  need  because  of  its  capability  as  a  flowsheeting 
package . 


In  Chapter  II  of  the  dissertation  a  description  of  the  program 
EROS  is  provided.   The  data  specification  and  modeling  considerations 
for  the  system  are  discussed  here  while  a  user's  guide  to  EROS  is 
presented  in  the  appendices.   The  derivation  of  a  solution  procedure 
is  illustrated  with  the  help  of  a  simple  example.   If  the  user  has  not 
provided  a  feasible  starting  point,  EROS  has  the  capability  of  dis- 
covering one.   The  algorithm  used  is  a  modified  version  of  that  pre- 
sented in  deBrosse  and  Westerberg  (1973).   The  optimization  strategy 
based  on  Westerberg  and  deBrosse  also  merits  attention  in  this  study 
because  of  its  usefulness  to  this  application.   The  different  applica- 
tions of  EROS  and  the  results  on  optimization  of  10  problems  of  varying 
sizes  are  also  shown. 

The  program  EROS  does  not  guarantee  a  global  optimum  because 
of  the  non-linear  and  multimodal  nature  of  the  problem.   In  Chapter 
III  an  algorithm  is  presented  to  establish,  often  quite  quickly, 
whether  the  discovered  optimum  is  global  or  not.   The  strategy  is 
based  on  finding  both  a  lower  (dual)  bound  for  the  optimum  of  the 
structured  system  and  an  upper  bound  on  this  lower  bound.   In  this 
chapter  the  effectiveness  of  the  algorithm  is  then  illustrated  by 
applying  it  to  three  small  heat-exchanger  network  problems. 

With  a  convenient  evaluation  and  optimization  package  such  as 
EROS,  it  is  very  tempting  to  use  it  to  perform  process  synthesis  via 
the  use  of  structural  parameters.   Using  this  approach  several 
alternate  configurations  are  imbedded  into  one  main  structure 
which,  on  optimization,  yields  the  required  test  configuration. 


Several  authors,  as  mentioned  earlier  in  this  section,  have  already 

experimented  with  this  idea.   However,  the  problem  to  be  solved  has 

to  be  formulated  very  carefully,  and  not  enough  attention  has  been 

paid  to  this  iact.   In  Chapter  IV  a  discussion  regarding  this  subject  is 

presented. 

The  conclusions  from  creating  a  system  such  as  EROS  and  recom- 
mendations for  further  investigations  form  Chapter  V  of  this  disserta- 
tion.  A  guide  to  the  use  of  the  program  is  given  in  the  appendices. 
The  aspects  of  EROS  presented  are  the  input  data  format,  the  inter- 
pretation of  the  output,  a  typical  computer  printout,  and  a  description 
of  all  the  subroutines. 


CHAPTER  II 
THE  SYSTEM  EROS 

II. 1.   Background 

In  order  to  evaluate  and  optimize  heat-exchanger  networks  it 
is  desirable  to  have  a  program  which  on  being  given  information  about 
the  configuration  and  stream  properties,  yields  all  the  required 
information  about  the  optimal  network.   Since  the  program  will  perform 
several  different  tasks,  it  would  be  very  attractive  and  in  many 
instances  necessary  for  i t  to  possess  the  following  features. 

a)  An  ability  to  set  up  solution  procedures. 

The  program  should  eliminate  or  reduce  com- 
putational recycles,  choose  the  decision  variables 
and  discover  the  order  for  calculating  the  various 
unknown  variables. 

b)  An  ability  to  obtain  a  feasible  starting  point. 

Often,  locating  a  feasible  starting  point  for 
optimization  is  not  a  simple  task.   In  order  to  save 
users  the  time  and  trouble  necessary  to  find  a  feasible 
starting  point,  the  program  should  be  capable  of  per- 
forming such  a  task  on  its  own. 

c)  Efficient  optimization  routines. 

These  would  be  required  for  selecting  the 
optimal  values  for  the  decision  variables  by 
searching  over  their  full  range. 

However,  the  optimization  of  a  system  such  as  this  one  raises 
certain  problems.   The  objective  function  is  highly  non-linear  and 
multimodal.   Also,  if  phase  changes  are  allowed,  continuous  derivatives 
cannot  be  obtained.   These  criteria  force  the  use  of  a  search  algorithm 
such  as  the  complex  method.   In  having  resigned  to  the  use  of  the 


complex  method  for  optimization,  further  demands  can  be  made  to  improve 
the  efficiency  of  the  approach  for  optimization.   If  in  the  process  of 
optimization  several  constraints  are  violated,  one  remedial  action  is 
the  use  of  penalty  functions,  but  this  modification  is  inefficient  and 
it  increases  the  number  of  iterations  required  for  convergence. 

Hence,  with  regard  to  optimization,  few  additional  features  would 
be  deemed  attractive. 

d)  The  ability  to  rederive  a  solution  procedure. 

When  a  constraint  is  violated,  the  program 
should  be  able  to  modify  the  equation  set  and  re- 
derive  a  solution  procedure  so  that  the  optimization 
may  be  continued. 

This  strategy  will  lessen  the  number  of  iterations  required 
for  convergence  as  compared  with  the  penalty  function  method.   However, 
it  is  essential  that  the  saving  in  computer  time  thus  incurred  com- 
pensates for  the  extra  time  required  in  rederiving  the  solution 
procedure . 

e)  The  use  of  restriction  as  a  solution  strategy. 

In  the  presence  of  a  large  number  of  decision 
variables,  some  of  them  are  set  to  zero  and  optimiza- 
tion is  performed  with  only  the  remaining  ones  as 
search  coordinates.   Optimization  is  considered 
complete  when  the  Kuhn-Tucker  (1951)  optimality 
conditions  are  satisfied  with  respect  to  the  decision 
variables  held  at  zero  value.   This  strategy  aids 
considerably  in  reducing  the  number  of  iterations 
during  optimization. 

EROS  is  the  author's  attempt  at  developing  an  optimization 
program  which  incorporates  all  the  features  discussed  above.   Figure  1 
illustrates  the  general  structure  of  the  EROS  system. 

The  approach  taken  is  to  model  each  unit  in  a  heat-exchanger 
network  functionally  by  writing  overall  material  and  energy  balances. 
The  unit  models  themselves  may  be  very  complex  internally,  but  the 


Objec tlve 
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Fig.  1.   Structure  of  an  Optimizing  Syster 


net  effect  at  the  flowsheet  level  is  that  each  unit  satisfies  overall 
heat  and  material  balances.   In  the  current  version  of  EROS  only  simple 
models  are  used.   Using  the  functional  equations,  the  modeling  of  such 
a  system  is  only  at  the  flowsheet  level,  and  a  solution  procedure  or 
the  order  in  which  these  equations  are  to  be  solved  is  found  along  with 
the  degree  of  freedom  to  be  chosen.   The  solution  procedure  sought  is 
one  that  will  eliminate,  if  possible,  computational  recycles  at  this 
level . 

The  above  approach  is  useful  because  the  equations  are  solved 
repeatedly  as  an  inner  loop  to  an  optimization  program.   As  illustrated 
in  Figure  1,  the  optimizer  directs  all  the  activity.   Its  primary  func- 
tion is  to  adjust  the  decision  variables  to  improve  the  objective 
function  0.   For  this  system  0  is  the  annualized  cost  of  the  equipment 
plus  the  annual  cost  of  buying  the  utilities  needed  such  as  steam  for 
the  purpose  of  heating. 

To  evaluate  0,  the  optimizer  supplies  the  block  labeled  "Solve 
Model  Equations"  with  the  values  it  wishes  to  try  for  the  decision 
variables  u.   The  remaining  problem  variables  x(u)  are  then  obtained 
by  solving  the  model  equations.   With  u  and  x(u)  values  available, 
constraint  violations  are  checked,  and  if  some  are  violated,  they  are 
identified  to  the  optimizer.   Assuming  none  are  violated,  the  units 
in  the  system  are  sized  and  an  annualized  cost  0  is  evaluated.   The 
optimizer  notes  this  cost  and  changes  the  decision  variable  values, 
with  the  aim  of  reducing  0.   This  calculation  sequence  is  repeated 
many  times  during  a  typical  optimization.   If  constraints  are  violated, 
special  action  is  taken,  which  tor  this  system  will  result  in  a  modi- 
fied set  of  model  equations  and  a  need  to  rederive  a  solution  procedure 


for  them  (this  approach  is  based  on  the  optimization  strategies  in 
deBrosse  and  Westerberg  (1973)  and  Westerberg  and  deBrosse  (1973)). 
The  modified  complex  optimization  algorithm  (Umeda  and  Ichikawa  (1971)) 
is  used  in  searching  for  the  decision  variable  values,  u. 

11. 2.  Data  Specification 

Adequate  data  must  be  supplied  to  the  computer  to  define  the 
problem.   A  problem  definition  will  require  the  following  input. 

1)  the  flowsheet  structure 

a)   connectivity  of  segments  and  units 

2)  the  unit  data 

a)  unit  type 

b)  any  equipment  parameter  specifications 

3)  the  stream  data 

a)  flow  and  temperature  specifications 

b)  physical  property  data 

c)  film  heat  transfer  coefficients 

d)  materials  specifications  and  other  cost 
parameters 

4)  the  segment  data 

a)  associated  stream  identifier  (i.e.,  what 
stream  this  segment  is  a  part  of) 

b)  any  specifications  imposed  on  flow  and 
temperature 

5)  general  user  specifications 

6)  guessed  set  of  inequality  constraints  to  be  held 

11 .3.  Modeling  Considerations 

The  modeling  of  a  network  will  be  illustrated  with  regard  to 
the  example  in  Figure  2. 

The  network  comprises  a  single  hot  process  stream  H  which  is 
split  and  used  to  heat  two  cold  process  streams  C   and  C„.   It  then 
merges  to  its  exit  conditions.   Streams  C   and  C„  are  heated  further 
by  steam  utility  streams  S   and  S,;.   The  network  has  four  heat 
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I  (11,1 
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I  -   500 
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p  2  3  5  6 
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For  Stream  III       Tja  -  T^  -  600°  neat  of  Vapomation  =•  800 
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vhere    Cost,    -  35   (/area.)    '    .    Cost,    and  Area. 
1  -  i  1  1 

axe   associated  u-lth  exchanger   1. 


Fig.   2.    An  Example  Problem 
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exchanger  (or  heater  and  cooler)  units,  2,  3,  5  and  6,  one  stream 
splitter  unit,  1,  and  one  mixing  unit,  4.   These  units  may  also  be 
referred  to  as  nodes.   All  the  heat  exchangers  are  assumed  counter- 
current.   The  streams  have  been  broken  up  into  segments,  of  which 
there  are  16  overall.   For  example,  stream  C1  enters  node  2  as 
segment  11.   It  exits  and  proceeds  to  unit  5  as  segment  12,  and 
finally  leaves  the  system  as  segment  13.   The  naming  scheme  should 
now  be  evident.   The  3  basic  units  used  are  shown  in  Figure  3.   The 
unit  models  are  written  functionally  by  writing  overall  material 
and  energy  balances. 


Unit 


F  =  aF  (I) 

F3  =  (1  -  a)F1                                (2) 

h2  =  hx  (3) 

h3  =  h1  (4) 


F4  «  F2  (5) 

(6) 
h2F2  +  hllFll  =  h4F4  +  h12F12  (?) 


(8) 

(9) 
h3F3  +  h14F14  =  h5F5  +  h15F15  (10) 


6     4b 

h,F,  =  h.F  +  h  F,  (12) 

6  6    4  4    5  5 
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cream  VyHcat-Exclirfiige 


d 


Fig.  3.   The  Network  Nodes 
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5  Fg  -  Fy  (13) 

F13  -  F12  (14) 

h7F7  +h12Fl2  =  h8F8+hl3F13  (15) 

(16) 

(17) 
h9F9+h15F15  =  h10F10  +  h16F16  (L8) 

In  addition  to  these  18  equations,  the  associated  inequality 
constraints  and  the  equipment  sizing  and  costing  relations  can  be 
written. 

The  basic  inequality  constraint  is  that  at  no  point  in  the 
exchanger  should  the  hot  stream  temperature  equal  or  fall  below  that 
of  the  cold  stream.   Referring  to  the  heat-exchanger  in  Figure  3a, 
this  constraint  is  usually  expressed  as 

T   >  T,  +  D,     where  D  is  minimum  allowable  approach 

1       4 

temperature . 
T2  >  T3  +  D 

However  the  temperatures  could  cross  over  internally  and  these 
constraints  may  not  be  adequate  to  detect  it,  particularly  when  a 
stream  passes  through  a  phase  change.   A  check  should  therefore  be 
made  at  several  points  along  the  exchanger  to  prevent  a  "crossover" 
of  temperatures. 

The  final  set  of  constraints  indicate  that  a  positive  heat 
transfer  must  occur,  that  is,  the  hot  stream  must  be  cooled  and  the 
cold  stream  heated.   These  are 
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T„  <  T.    and    T-  <  T.  . 
2     1  3     4 

All  Che  constraints  associated  with  a  typical  exchanger  such 
as  the  one  shown  in  Figure  4  are  listed  in  Table  1  with  their  ap- 
propriate code  so  that  they  can  be  precisely  identified  by  numbers. 
No  constraints  are  written  for  the  splitter  and  mixer  units. 

The  sizing  calculation  for  an  exchanger  is  to  evaluate  its 
area.   This  calculation  can  be  very  involved,  but  for  design  purposes 
may  in  fact  be  simplified  by  assuming  film  coefficients  based  on  the 
fluid,  whether  it  is  heating  or  cooling  and  whether  it  is  boiling  or 
condensing  (see  Perry  and  Chilton  (1973)).   The  exchanger  may,  again, 
for  simplified  design  purposes,  be  considered  to  operate  in  zones  as 
indicated  in  Figure  4.   Each  zone  is  then  sized  and  costed  using  the 
appropriate  film  coefficients  and  temperature  driving  forces. 

To  place  a  crude  cost  on  the  exchanger  one  might  use  an  equation 
of  the  form  (see  Guthrie  (1969)). 

COSt  =  fMfP(aAm) 

Where  f   is  a  materials  factor,  fD  a  pressure  factor  and  A  the 
M  ' 

area  of  a  zone  within  an  exchanger.   The  zones  are  sized  and  costed 
differently  because  of  the  different  types  of  heat-exchange  duty. 
The  terms  a  and  m  are  constants,  with  m  being  about  0.6  to  0.8. 
Constant  costs  are  assumed  for  the  splitter  and  the  mixer  units. 

The  last  source  of  equations  is  the  evaluation  of  physical  prop- 
erties.  The  system  must  be  able  to  convert  from  stream  temperature  (and 
vapor  fraction  for  a  pure  component  in  two  phase  region)  to  enthalpy 
and  vice  versa.   A  cooling  curve  may  be  provided  for  the  stream 
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TABLE    1 
CONSTRAINT   REPRESENTATION   FOR  AN   EXAMPLE   NODE    3    IN   FIGURE   4 


Exterior 

Cggggent 

Constraints 

Code 

approach 

T   s  T     +   D 

(NODE 

*  10)  +  i  =  31 

approach 

Tjit  +  D 

(NODE 

*  10)   +   2  =  32 

VT2 

(NODE 

•  10)  +  3  *  33 

T^T3 

(NODE 

*  10)  +  U  =  34 

F1IR  =  lower  bound 
LUi       for  How  F 

F12FLB- 

(NODE 

*  10)  +  5  =  35 

F.  _.  =  upper  bound 
1U^       for  flow   F, 

F     2  F 

OB       1 

(MODE 

»  10)  +  6  =  36 

?-.,  u  -  lower  bound 
^         for  flow  F 

F   2  F 
r3       3  LB 

(NODE 

*  10)  +  7  =  37 

F3DB  =  upper  bound 
J           for  flow  r 

F       a  F 
3  OB       3 

Interior 

(NODE 

»  10)  -t    8  =  38 

Type 

Constraints 

Representation 

approach 

VW  d 

(NODE 

*  1000)  +  1  =  3001 

approach 

WT6  +   D 

(NODE 

*  1000)  +  2  =  3002 

approach 

.  VW  ° 

(NODE 

*  1000)  +  3  =  3003 

approach 

W° 

(NODE 

*  1000)   +  A  =  3004 

In  addition  there   could 

be   constraints  asr 

jciated  with   any   stream  seement. 

Example:     Segaent  3,   No 

le  3 

Coaiuent 

Constraints 

Code 

H,.  R  -   lower  bound 
J            for  enthalpy  h\ 

K3iH3LB 

-   (SEG  « 
+   1) 

1000  +    NODE   »  10 
^-3031 

H^rm  =  uPPer  bound 

^OB*^ 

-    (SEG   " 

4.     2\ 

1OC0  +   NODE  *  10 

for  enthalpy  H 
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if  the  stream  is  assumed  to  be  at  a  constant  pressure.   Figure  5 
illustrates  a  cooling,  curve,  where  T  and  T   are  the  dew  and 
bubble  point  temperatures  respectively. 

Properties  such  as  thermal  conductivities  and  densities 
should  also  be  provided  if  the  film  coefficients  are  to  be  deter- 
mined from  correlations.   If  for  design  purposes  typical  values  for 
film  coefficients  are  to  be  used,  these  properties  will  not  be  needed. 

II  .4  .   Deriving  ^-h£  Solution  Procedure  and_ 
Solving  Model  Equations 

Consideration  will  now  be  given  to  developing  a  solution 
procedure  and  then  solving  the  example  problem.   First,  the  system 
must  gather  together  the  necessary  equations,  or  at  least  establish 
their  structure,  so  that  a  solution  procedure  may  be  prepared.   The 
desired  solution  procedure  should  eliminate  all  recycle  loops  in  the 
computations  if  possible  or  minimize  their  number. 

The  initial  solution  procedure  will  ignore  all  but  the  user 
specified  inequality  constraints.   Thus  the  system  sets  up  the  18 
heat  and  material  relations  shown  in  the  last  section.   Assuming  that 
the  user  has  requested  that  constraints  55  and  33  be  included,  the 
additional  relations  required  are 


F7  =  F7LB  +  °55 


T3  -  T5  +  a33 


(19) 

(20) 


where  F     is  the  lower  bound  on  F7. 
/LB  / 

The  inequality  constraints  have  been  converted  to  equality  by 
the  slack  variables  o  ,  and  Q„„  which  are  required  to  be  nonnegative. 


18 


T      "   Bubble   point    tecrpcratur 
T     -  Dew  point    temper  a  Cure 


Fig.    5.      A  Typical   Cooling   Curve 


19 


When  the  solution  procedure  is  derived  o   and  a   will  be  required 

to  be  decision  variables  with  an  initial  value  of  zero.   In  this  way 

F-,  and  T„  will  be  forced  to  equal  to  F    and  T   respectively. 
/       3  /Lis       j 

The  relation  (20)  is  in  terms  of  temperatures  rather  than 
enthalpies.   Hence  the  following  relationships  should  also  be  added 

T3  -  f(h3)  (21) 

T5  =  f(h5)  (22) 

The  system  can  implicitly  account  for  equations  (20) ,  (21) 
and  (22)  by  the  single  expression 

h3  =  f(h5,o33)   .  (23) 

The  incidence  matrix  can  now  be  created.   However  its  size  can 
be  significantly  reduced.   Note  that  a  large  number  of  equations  simply 
equate  one  variable  to  another.   Equations  (3),  (4),  (5),  (6),  (8),  (9). 
(13),  (14),  (16)  and  (17)  are  precisely  of  this  form.   These  equations 
will  automatically  be  satisfied  if  the  values  of  variables  so  equated 
are  stored  in  a  common  storage  location.   Hence  these  equations  may 
be  deleted  and  the  two  variables  occurring  in  each  one  of  them  may  be 
merged  to  a  single  one. 

Many  of  the  variables  in  the  incidence  matrix  are  in  fact 
specified.  Let  the  following  be  specified  in  data  input  for  the 
example  in  Figure  2 

Flows       Fi'Fii'Fi4'F7LB 

Enthalpies   ^V'? 'h8'Vh10'hll 'h13>h14'h16 
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These  specified  variables  along  with  the  slack  variables  o   and  <J 
(required  to  be  decision  variables)  may  be  eliminated  in  the  incidence 
matrix.   The  resulting  and  much  reduced  matrix  is  illustrated  in 
Table  2. 

A  modification  of  Lee,  Christensen  and  Rudd  (1966)  algorithm 
is  applied  to  determine  the  solution  procedure.   The  solution  proce- 
dure that  results  on  the  application  of  this  algorithm  to  the  incidence 
matrix  of  Table  2  is  shown  in  Table  3.   The  variables  listed  are  cal- 
culated from  the  corresponding  equations  in  the  order  indicated.   Note 
that  there  is  an  iteration  loop  involving  the  single  'tear'  variable 
F   (from  steps  A  to  9).   F_(=F.)  appears  in  equation  (12)  and  the 
iterations  between  steps  4  and  9  are  continued  until  the  value  of  F.; 
guessed  at  step  4  is  essentially  the  same  as  the  value  of  F^  calculated 
from  equation  (12)  in  step  9. 

The  execution  of  the  solution  procedure,  that  is,  calculation 
of  the  variables  from  the  equations  assigned  to  them,  is  termed 
"Solve  Model  Equations"  in  Figure  1.   Corresponding  to  every  unit,  a 
subroutine  is  required  to  calculate  any  variable  involved  in  the  heat 
and  mass  balances  of  the  particular  unit.   These  subroutines  may  be 
supplied  by  the  user  for  more  sophisticated  models. 

II.  5.   Starting  the  Problem:   Finding  a  Feasible  Point 

If  the  user  has  not  provided  any  information  to  aid  in  obtain- 
ing a  feasible  starting  point,  a  modified  version  of  an  algorithm  by 
deBrosse  and  Westerberg  (1973)  is  used.   As  mentioned  earlier,  a  signif- 
icant effort  would  be  required  on  the  part  of  a  user  to  provide  a 
feasible  starting  point  if  computational  loops  are  involved  in  solving 
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TABLE  2 
THE  INCIDENCE  MATRIX 


Variables 

Equations 

tx 

;- 

(1) 

(2) 

(7) 

(10) 

(11) 

(12) 

(15) 

X 

(18) 

(19) 

X 

(23) 

-<f 

IT\ 

c\j 

LPv 

t  1 

r-i 

IT\ 

xf1 

-*0 

a 

X 

X 

X 

X 

X  XX 

XX  x 

XX  XX  X 
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TABLE  3 
SOLUTION   PROCEDURE  FOR  TOE  PROBLEM   IN   FIGURE  2 


Decision   Variables 
Variable 


i.  y=v 


n  ,  _  .  it 


2-    ^ 

3.     h5 


4.      Guess   Fp(=F.  ) 


5.     a 


6.      F   (=F   ' 


7.      h, 


8.      F, 


'•  *-2<=y 


10.      H 
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ii.   f9(=f10: 


55'      33 
Equation 

(19) 

(15) 
(23) 

(1) 
(2) 


HI) 

(12) 
(10) 
(18) 
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model  equations.   Computational  loops  are  almost  inevitable  in  complex 
networks.   However,  the  user  does  have  the  option  of  providing  a  feasible 
starting  point. 

The  deBrosse  and  Westerberg  (1973)  algorithm  uses  an  indirect 
approach.   It  hypothesizes  that  a  subset  of  constraints  has  no  feasible 
region  and  then  attempts  to  verify  the  conjecture.   If  successful,  the 
subset  is  identified  as  infeasible  and  obviously  no  feasible  point 
exists.   If  unsuccessful,  either  a  new  hypothesis  can  be  generated  or 
the  algorithm  has  indirectly  found  a  feasible  point. 

The  feasible  pq i n t __al go r i t hm 

In  order  to  demonstrate  the  feasible  point  algorithm,  the 

following  definitions  are  necessary. 

E   =  The  equation  set  without  any  inequality  constraints  present 
as  equalities  (heat  and  material  balances  only)  . 

E   =   The  current  equation  set  (E.,  +  inequality  constraint  (s) 
present  as  equalities  with  slack  variables) . 

C   =  The  set  of  all  newly  violated  constraints. 

C.   =  The  violated  constraint  encountered  first. 

V  =  The  current  set  of  inequality  constraints  in  the  equation 
C     set:   (Ec  -  EQ). 

H   =  A  set  of  (m  +  1)  constraints. 

V  ,V  , . . .V   .   =   subsets  formed  by  removing  one  constraint  at 

a  time  from  H  (for  example,  V-i  is  obtained  by 
removing  the  first  constraint  in  H) . 

'0'  =   null  set. 

Degree  of  freedom  (as  defined  for  this  work)  =  The  number  of  original 
variables  (excluding  slack  variables)  in  problem  less  number  of  con- 
straints in  E  . 
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Procedure  'A'  represents  the  following  sequence  of  operations. 

1)  Set  all  the  decision  variables  to  a  value  of 
zero . 

2)  Solve  the  model  equations. 

3)  Check  for  constraint  violations. 


Step  1      Set  E   =  E„,  V  =  '0' 
c     0    c 

Step  2      Determine  the  solution  procedure  based  on  struc- 
tural considerations  of  the  equation  set.   If  a 
solution  procedure  can  be  determined,  go  to  Step  4. 

Step  3      Attempt  unsuccessful.   Exit. 

Step  4      Execute  Procedure  'A' .   If  any  constraints  are 
violated,  go  to  Step  6. 

Step  5      Feasible  starting  point.   Exit. 

Step  6      Set  E   =  [E  ,C,  ] 
c      c   1 

Step  7      Determine  the  solution  procedure.   If  a  solution 
procedure  cannot  be  found,  go  to  Step  11. 

Step  8      Execute  Procedure  'A'.   If  the  problem  does  not 
prove  solvable,  go  to  Step  11.   If  no  constraints 
are  violated,  go  to  Step  5. 

Step  9      If  the  degree  of  freedom  is  not  zero,  go  to 
Step  6. 

Step  10     Set  H  =  [V  ,C  ],  i  =  1.   Co  to  Step  12. 

Step  11     Set  i  =  1,  H  =  V 

Step  12     Set  E   =  [E„,V.],  V.  +   V   from  Step  10. 
c      0   l     l     c 

Step  13     Determine  the  solution  procedure.   If  it  cannot 
be  determined,  go  to  Step  17. 

Step  14      Execute  Procedure  'A' .   If  the  problem  cannot  be 

solved  go  to  Step  17.   If  no  constraints  are  violated, 
go  to  Step  5. 

Step  15     If  C  Pi  H  f    '0',  go  to  Step  17. 

Step  16      If  the  degree  of  freedom  =  0,  go  to  Step  10.   If 
not,  go  to  Step  6. 
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Step  17  Set  i  -  i  +  I.  If  i  >  m  +  1,  go  to  Step  3  (the 
set  of  constraints  in  H  do  not  allow  for  a  feasible 
region).   If  i  <•   m  +  1,  go  to  Step  12. 


The  main  departures  from  the  deBrosse  and  Westerberg  (1973) 
algorithm  are  as  follows: 

1)  In  this  work,  the  set  H  is  restricted  to  containing 
less  than  or  equal  to  (n  +  1)  constraints  where  n 
represents  the  number  of  decision  variables. 

2)  Both  the  sets  V   and  H  are  created  in  a  different 
manner.   In  this  work,  one  constraint  is  added  at  a 
time  to  build  up  the  set  V  ,  and  H  is  created  either 
when  there  are  already  n  constraints  in  V   and  an 
additional  constraint  is  violated,  or  when  a  struc- 
tural inf easibility  occurs. 

3)  "No  Solution"  options  (in  Steps  8  and  14  where  the 
problem  is  not  solvable)  are  not  treated  in  as 
rigorous  a  fashion  as  in  deBrosse  and  Westerberg 
(1973) .   The  method  here  is  the  same  unless  several 
problems  corresponding  to  the  same  "hypothesis" 

set  H  lead  to  no  solution  in  Step  14.   The  algorithm 
here  may  terminate  unsuccessfully  at  Step  17  whereas 
that  of  deBrosse  and  Westerberg  might  continue  with 
alternate  and  reduced  sets  H.   This  option  is  quite 
complex  and  was  never  found  necessary  in  EROS. 
Hence  it  was  never  included. 

The  algorithm  can  now  be  applied  to  the  example  shown  in  Figure 

2.   (Note  that  this  problem  involves  8  equations  in  10  unknowns;  thus 

two  degrees  of  freedom  exist  at  the  start.) 


EQ  =  [(1),  (2),  (7),  (10),  (11),  (12),  (15),  (18)] 

(all  the  equations  from  Table  2  except  the  last  two, 
(19)  and  (23)) 

1)  E   =  E..,  V   =  '0* 

c     0    c 

2)  Solution  procedure  derived 

4)   Procedure  'A'  executed 

Constraint  Violation  =  33. 
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E   =  [En,  33] 
c      U 

Solution  procedure  derived 

Procedure  'A'  executed 

Constraint  Violation  =  25 

Degree  of  freedom  =  1 

Ec  =  [EQ5  33,  25] 

Solution  procedure  derived 

Procedure  'A'  executed 

Constraint  Violation  =  22 

Degree  of  freedom  =  0 

H  =  [33,  25,  22],  V  =  [33,  25],  V2  =  [22,  33],  V3  =  [22,  25; 

Ec  =  [EQ,  V2]  -  [EQ,  22,  33] 

Note:   V   =  [33,  25]  is  V   in  the  previous  iteration  and 
it  need  not  be  considered  again 

Solution  procedure  derived 

Procedure  'A'  executed.   Constraint  Violation  =  55 

C  0  H  =  '0' 

Degree  of  freedom  =  0 

H  =  [22,  33,  55],  V   =  [22,  33],  V2  =  [55,  22],  V3  =  [55,  33 j 

Ec  =  [EQ,  V2]  =  [EQ,  55,  22] 

Note:   V  =  [22,  33]  is  V   in  the  previous  iteration  and 
is  eliminated  from  consideration 

13)  Solution  procedure  derived 

14)  Procedure  'A'  executed 
No  constraint  violation 

5)   Feasible  starting  point.   Exit. 
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II . 6.   The  Optimization  Strategy 

The  optimization  strategy  is  modeled  after  the  algorithm 
presented  in  Westerberg  and  deBrosse  (1973)  .   The  algorithm  is  in- 
voked once  a  feasible  point  is  available. 

The  sets  of  inequality  constraints  are  divided  into  three  sets. 

V    =    The  set  of  inequality  constraints  being  held  as 
equality  constraints . 


V  =   The  set  of  constraints  present  in  the  equation 

set  as  equality  constraints  with  the  difference 
that  their  slack  variables  are  used  as  search 
coordinates. 

V  =   The  set  of  all  remaining  constraints. 

V  =    [Vm,V„]   (Note:   V   is  as  defined  in  the  previous 
c        T   R  c 

section. ) 

Solution  procedures  are  modified  as  inequality  constraints  are 
moved  from  one  set  to  another.   Adding  constraints  in  the  set  being 
held  tends  to  aid  the  optimization  process  by  reducing  the  dimension 
of  search  space  for  what  is  usually  a  marginal  or  no  added  burden  in 
solving  an  enlarged  set  of  equations. 

As  the  optimization  proceeds,  the  values  of  all  variables,  in- 
cluding all  slack  variables,  are  stored  for  the  point  that  yields  the 
best  value  for  the  objective  function.   Hence,  even  when  V^  is  changed, 
optimization  can  be  and  is  started  at  the  best  point  discovered  up  to 
that  moment.   This  modification  makes  a  significant  improvement  to 
the  Westerberg  and  deBrosse  (1973)  method.   Figure  6  illustrates  the 
typical  dilemma  faced  by  the  Westerberg  and  deBrosse  optimization 
algorithm  when  stepping  from  a  current  best  point,  point  'e',  through 
one  or  more  inequality  constraints  to  point  'f'.   At  point  'f'  the 


:-h 


Fig.  6.   Keeping  Track  of  the  Best  Point  'e' 
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constraint  g1  is  detected  as  being  violated.   The  algorithm  will 
respond  by  changing  the  solution  procedure  so  that  the  slack 
variable  0   for  g   becomes  a  decision  variable.   The  other  decision 
variable  will  be  either  x   or  x  .   There  are  several  options  now  as 
to  where  the  optimization  may  be  started.   The  algorithm  could  hold 
x.  or  x„  (whichever  is  selected  as  the  decision  variable)  at  its 
current  value  and  find  the  point  where  a      is  zero,  leading  to  point 
P,  or  P„  respectively.   Alternatively  it  could  attempt  to  locate  P„ 
by  searching  along  the  direction  leading  from  'e'  to  'f'.   All  of 
these  options  can,  and  often  do,  lead  to  a  next  point  which  has  a 
higher  and  thus  worse  value  for  the  objective  function.   By  saving 
all  the  variable  values  for  the  best  point,  the  search  can  always 
start,  even  after  developing  a  new  solution  procedure,  from  that 
point,  that  is  from  point  'e'.   This  change  reduces  cycling  because 
a  change  in  the  solution  procedure  cannot  lead  to  a  point  that  is 


The  actual  search  strategy  used  is  the  modified  complex 
method  (Umeda  and  Ichikawa  (1971)).   The  complex  method  is  consid- 
ered suitable  because  gradients  are  not  required.   The  treatment 
of  phase  changes  creates  discontinuities  in  first  derivatives  of 
functions . 

The  optimization  Algorithm 

The  notation  followed  is  the  same  as  that  in  the  last  section. 

Figure  7  gives  the  structure  of  this  algorithm. 

Step  1     Obtain  a  feasible  starting  point. 

Step  2     Set  VT  =  V  ,  VR  =  '0'.   Apply  the  Kuhn-Tucker  (1951) 
conditions  as  follows.   Perturb  each  slack  variable 
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Unsuccessful    Exit 


A   -   Horro/ll    mode    loop 

B   -   Find    scorch   coordin.ntrs,    nbnonal    node    loop 

C  -   fVCul  n    to   ii.-imal    n  iitc 


Fig.    7.      Structure   of   the   Optimization  Algorithm 
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corresponding  to  constraints  in  V^,  one  at  a 
time,  away  from  zero.   If  an  improvement  in 
the  objective  function  results,  move  the  cor- 
responding constraint  from  V   to  V  . 

A.  If  any  solution  procedure  fails  numerically  during 
this  step,  go  to  Step  11. 

B.  If  on  completion  of  this  step,  the  number  of  con- 
straints in  Vj  is  equal  to  the  number  of  constraints 
in  Vc,  and  the  degree  of  freedom  is  zero,  go  to  Step 
18. 

C.  Otherwise  continue. 


The  optimization  is  back  in  its  normal  mode  if  return  is  to  Step  3 
or  Step  7. 

Step  3        Perform  an  optimization  (using  the  complex  algo- 
rithm) .   Use  as  search  variables  the  slack  variables 
corresponding  to  inequality  constraints  in  VR  along 
with  any  additional  problem  variables  still  retained 
as  decision  variables.   Search  only  over  nonnegative 
values  of  the  slack  variables. 

A.  If  at  any  point  the  solution  procedure  fails  numeri- 
cally, go  to  Step  11. 

B.  Otherwise,  when  the  optimization  algorithm  exits 
normally,  continue. 

Step  4      An  optimization  has  just  been  completed. 

A.  If  one  or  more  constraints  are  violated,  go  to 
Step  7. 

B.  Otherwise  continue. 


Step  5        Set  V™  ,  ,  =  V,,,.   Apply  the  Kuhn-Tucker  conditions 
as  done  in  Step  2. 

A.  If  any  solution  procedure  fails  numerically,  go  to 
Step  11. 

B.  Otherwise  on  completion  continue. 


m 


Step  6  A.    If  V_  =  V,T,  .,,  go  to  Step  18. 
1     Iold 

B.    Otherwise  repeat  from  Step  3. 

One  or  more  constraints  not  in  the  current  set  V   has  been  violated. 

c 

Step  7  A.    If  the  degree  of  freedom  is  not  zero,  set 
E  =  [E  ,C ,]  and  go  to  Step  9. 

B.  Otherwise  find  the  set  of  constraints  R  from  V 
such  that  each  constraint  in  R  could,  from  struc- 
tural considerations  of  the  equations,  be  traded 
for  C  .   If  R  =  '0',  go  to  Step  10. 

C.  Find  the  constraint  C,  in  set  R  having  the 
largest  value  for  its  corresponding  slack 
variable.   If  that  value  is  zero,  go  to 
Step  10. 

D.  Otherwise  continue. 

Step  8      Replace  C  with  C   in  V  . 

Step  9        For  any  slack  variable  corresponding  to  a  constraint 
in  V^  and  having  a  value  of  zero  at  the  current  best 
point,  transfer  the  corresponding  constraint  from  Vp 
to  VT.   Set  V   =  [VD)Vn>],  E  =  [E„,V  ]  and  develop 

1  .  .     c       K  ,  *      C       UC 

a  new  solution  procedure. 

A.  If  a  solution  procedure  cannot  be  developed,  go 
to  Step  11. 

B.  Otherwise  return  to  Step  3. 


The  degree  of  freedom  is  zero  and  either  (a)  the  values  of  all  the 

slack  variables  in  V..  are  zero  or  (b)  the  set  R  is  empty. 
K 

Step  10      Set  H  =  [V  ,C  ],  set  i  =  1,  and  go  to  Step  12. 

The  current  set  of  equations  are  (or  appear  to  be)  numerically 
inconsistent. 

Step  11     Set  H  =  [V  ] ,  set  i  =  1. 
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Step  12     Set  E   =  [E„.V.],  V.  +   V   (this  test  is  only  relevant 
cOiiC 
if  entry  was  originally  from  Step  10  above)  .   Set  VR 

equal  to  the  set  of  all  constraints  in  V.  whose  slack 

variable  values  are  greater  than  zero.   Set  V.,,  equal 

to  all  remaining  constraints  in  V,. 

A.  If  VR  =  '0'  ,  go  to  Step  13. 

B.  Otherwise,  determine  a  new  solution  procedure. 

1.  If  one  can  be  found,  go  to  Step  14. 

2.  Otherwise,  go  to  Step  17. 

If  entry  to  Step  12  was  via  Step  10  originally,  a  degenerate  problem 

is  at  hand  as  V  and  C,  (NT7  +  1  constraints)  are  all  satisfied  in  a 
c      1   Vc 

subspace  of  dimension  NTI  . 
Vc 

Step  13     Determine  a  solution  procedure. 

A.  If  one  cannot  be  determined,  go  to  Step  17. 

B.  Otherwise,  perturb,  one  at  a  time,  the  slack 
variables  corresponding  to  the  constraints  in 

V 

1.  If,  on  any  perturbation,  constraints  are 
violated  such  that  HOC  i    '0',  go  to  Step  17. 

2.  If,  on  any  perturbation  of  a  slack  variable, 

an  improvement  is  obtained  in  the  objective 

function,  put  the  corresponding  constraint 

into  V„ . 
K 

C.  If,  on  completion  of  the  perturbations  in  Step  B, 
the  set  V_  =  '0',  go  to  Step  18. 

D.  Otherwise  continue. 

Step  14       Perform  an  optimization.   Use  as  search  variables 
the  slack  variables  corresponding  to  inequality  con- 
straints in  VR  along  with  any  additional  problem 
variables  still  retained  as  decision  variables.   Use 
the  solution  procedure  derived  in  Step  13.   Search 
only  over  nonnegative  values  of  the  slack  variables. 
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A.  If  at  any  point  the  solution  procedure  fails 
numerically,  go  to  Step  17. 

B.  When  the  optimization  algorithm  exits  normally 
and  if  a  constraint  violation  occurs,  then 

1.  If  HflC  /  '0',  go  to  Step  17. 

2.  Otherwise,  HOC  =  '0'.   Go  to  Step  7. 

C.  Otherwise  when  the  optimization  algorithm  exits 
normally  and  no  constraints  are  violated,  continue. 

Step  15       Set  V,„  ,  ,  =  V™.   Apply  the  Kuhn-Tucker  conditions 
,  fo]d     T     ' 

(as  done  xn  Step  2)  . 

A.  If  any  solution  procedure  fails  numerically,  go 
to  Step  17. 

B.  Otherwise  continue. 

Step  16A.    If  VT  =  VTold,  go  to  Step  18. 
B.    Otherwise,  return  to  Step  14. 

■k 

The  current  solution  procedure  failed  or  led  to  an  immediate 
constraint  violation  of  the  constraint  dropped  in  set  H. 

Step  17     Set  i  =  i  +  1. 

A.  If  i  >  m  +  1,  go  to  Step  19. 

B.  Otherwise  repeat  from  Step  12. 

Normal  exit.   Optimization  attempt  apparently  successful. 
Step  18     Optimization  complete.   Exit. 

Abnormal  exit.   Optimization  attempt  aborted. 

Step  19     Optimization  failed.   Exit. 

Application  of  the  strategy  to  the  Example  Problem  in 
Figure  2. 
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1)  Feasible  starting  point.   V  =  [55,22]   (from  the 
last  section) . 

2)  VR  =  [22] 

3  and  4)  Optimization  with  the  complex  method. 

7)  Constraint  violated  =  35,  R  =  [22],  and  oJ0    >    0. 

8)  Replace  constraint  22  with  constraint  35. 

9)  VR  =  [35],  VT  =  [55].   Ec  =  [EQ,VR,VT]. 

3)  Optimization  with  the  complex  method. 

4)  Optimization  complete. 

5  and  6)   Vm  1 ,  =  [55].   V^  =  [55]  after  Kuhn-Tucker  test. 
iold  T 

18)   Optimization  successful.   Exit. 
At  the  optimum,  the  following  values  resulted: 

F   =  F.  =  7350,  F„  =  Fc  =  2650,  T, „  =  500°,  T, c  =  300°, 
2     4  3     5  12  15 

F7  =  0,  Fn  =  2500,  T.  =  392°,  T   =  423°  and  0  =  3538 

/     y       4       5 

II. 7.   Special  Uses  of  EROS 

Using  existing  Heat-Exchangers.   The  program  may  be  used  to 
analyze  a  network  where  some  of  the  exchangers  are  already  specified. 
A  special  effort  must  be  made,  however,  to  make  use  of  these  exchangers. 
It  is  assumed  that  these  exchangers  are  available  at  no  cost.   In 
Figure  8,  an  exchanger  with  an  area  of  A  has  been  specified.   On 
analysis,  however,  it  is  discovered  that  an  exchanger  with  an  area  A 
is  required  at  that  particular  site  in  the  network.   In  the  program, 
the  costs  assumed  for  different  conditions  of  A   and  A   are  shown  in 
Figure  8.   The  physical  significance  of  1)  is  that  a  by-pass  will  be 
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Specified 
Exchanger 

Require 
Exchang 

1) 

2) 

If  Aj«  A, 

If  A,  >  A, 

T03t 
Cost 

Cost  associated 
with  an  exchange 
having    the    axes 


Fig.    8.      Using   an  Existing  Exchanger 
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used  in  the  specified  exchanger.   2)  implies  that  an  exchanger  with 
area  =  A„  -  A  must  be  purchased  in  addition  to  the  exchanger  already 
present . 

Reliability  analysis.   The  program  has  been  extended  to  permit 
its  use  in  quick  reliability  studies.   The  reliability  studies  will 
be  demonstrated  with  the  help  of  an  example.   Figure  9a  represents 
a  network  that  is  operational  under  normal  conditions.   Cold  process 
streams  C, ,  C  ,  and  C   are  heated  to  their  final  temperatures  with 
the  help  of  a  hot  process  stream  H,,  and  a  flue  gas  H  .   The  flow 
rate  and  the  outlet  temperature  of  stream  H.  are  undefined  but  are 
required  to  be  within  specified  bounds.   Assume  that  two  abnormal 
occurrences  take  place  separately,  for  certain  periods  of  the  year, 
namely , 

1)  Stream  H„  is  unavailable. 

2)  Heating  of  stream  C   is  no  longer  required. 

The  aim  now  is  to  find  an  optimal  network  such  as  the  one 
shown  in  Figure  9a,  fully  provided  for  to  meet  the  contingencies  with 
the  aid  of  by-passes  in  the  exchangers  or  with  the  aid  of  auxiliary 
exchangers.   The  designer  is  permitted  to  allow  a  change  in  the 
flow  of  stream  H  and  its  outlet  temperature  provided  they  stay  within 
their  specified  bounds.   In  the  case  of  failure  mode  2),  a  cooling 
utility  stream  may  be  used  to  cool  stream  C„  to  230°  so  that  it  may 
be  recycled  to  be  heated  to  250°  if  it  is  desired  to  maintain  a 
semblance  to  the  normal  operation.   (Flow  rate  of  C„  can  be  allowed 
to  vary  between  0  and  70,000  in  this  instance.) 
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F    -    '•0,000 
T    -    150° 


F    »    20,000 
T    -    300° 


1      > 
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F    -    70,000 
I   -    30° 
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(c) 

*  Uaspeclfled 

Fig.    9.      Reliability  Analysis 
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Figure  9b  represents  the  network  when  stream  H   is  unavailable, 
and,  Figure  9c  when  stream  C„  is  not  required  to  be  heated.   Additional 
user  specifications  to  those  shown  in  Figure  9  are  presented  in  Table  4. 
In  order  to  find  the  optimal  operating  system,  networks  in  Figure  9a,  9b 
and  9c  are  optimized  together.   The  objective  function  0  is  given  by 

0  =  ch1fh1  +  ch;fh;  +  ch-'fh^  +  cuc"Fuc"  + 

5 

I    (cost  of  exchanger  area  at  site  i) 
i=l 

where  C.  and  F.  represent  the  cost  coefficient  and  the  flow  rates  of 
l      l 

stream  i,  respectively.   The  cost  coefficient  C.  should  reflect  the 

l 

expected  fraction  of  the  year  that  the  network  is  in  the  particular 
state  being  represented.   For  example,  for  the  problem  in  Figure  9 
it  is  assumed  that  the  networks  in  Figure  9a,  9b  and  9c  are  opera- 
tional 77%,  11.5%,  and  11.5%  of  the  time  in  a  year,  respectively. 

Hence  if  C    is  0 . 1 ,  then  C  ,  and  C  „  are  0.015  each. 
1  1       1 

The  cost  of  exchanger  area  at  a  site  will  be  illustrated  with 

the  help  of  an  example.   At  the  site  2,  the  exchangers  of  different 

sizes  required  in  Figure  9  are  A  ,  A  ,  and  A  ,,.   Assume  that  A  ,  is 

smallest  and  A  ,,  the  largest  area. 

The  cost  of  exchanger  area  at  site  2  is  defined  as 

35[(A2,)U'6  +  (A2  -  A2,)°-6  +  (A2„  -  A2)°-6] 

This  manner  of  costing  areas  in  doing  reliability  analysis 
appears  to  be  a  good  formulation  of  the  real  system.   If,  because  of 
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TABLE  4 
ADDITIONAL  USER  SPECIFICATION  FOR  THE  PKOBEEM  IM  FIGURE  9 


Stream 


UC" 


Flow  Rate 


Lower  Bound  on  Flow 


50,000  50,000        50,000  0 


Upper  Bound   on  Flow  200,000       200,000     200,000     100,000 


Inlet  Temp 


600 c 


600°  600°  100° 


Outlet  Temp 


150° 


Lower  Bound   on  Outlet   Temp       190°  100°  190°  150C 


Upper  Bound  on  Outlet  Temp       600°  600°  600°  150c 


Cost  Coefficient 


0.10  0.015       0.015  0.008 


Unspecified. 

Let  U.  represent  the  heat  transfer  coefficient  for  exchanger  I 

ui ;  V  =  7CJ0 

U2   =    V    "  l,2"   "  Ul1  •2l 

U3  =    Uy  U,»  -    562.5 

U,   -•    U    '    =  U  "   --,   700 

U  •'  =  225 
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some  departure  from  the  normal  mode,  more  exchanger  area  is  required 
at  a  particular  site,  then  one  must  pay  for  the  auxiliary  exchanger. 
If  the  exchanger  area  required  is  more  for  the  normal  mode,  then  an 
auxiliary  exchanger  is  already  in  effect.   The  results  obtained  on 
optimization  of  the  system  in  Figure  9  are  shown  in  Table  5. 

II. 8.   Results  and  Discuss ion 

A  typical  application  of  EROS  to  a  heat  exchanger  network  is 
illustrated  by  the  example  in  Figure  10,  the  stream  specifications 
for  which  are  shown  in  Table  6.   Stream  I  is  a  flue  gas  and  streams 
IV  and  VIII  are  utility  streams.   The  flow  rates  for  these  streams 
are  not  defined  but  are  required  to  be  within  specified  bounds. 
Some  of  the  streams  in  the  network  are  also  multicomponent  and  are 
characterized  by  a  dew  and  a  bubble  point.   Thus  the  network  is 
analyzed  to  ensure  that  minimum  approach  temperature  violation  does 
not  occur  inside  the  exchanger  owing  to  discontinuities  at  the  dew 
and  bubble  points.   In  fact,  at  the  optimal  solution  for  the  example 
in  Figure  10,  the  minimum  approach  temperature  constraint  at  the  bubble 
point  of  stream  III  is  operative  inside  exchanger  5.   The  program  can 
also  evaluate  variables  that  already  exist,  fully  or  in  part.   For 
example,  exchanger  3  was  assumed  present  with  an  area  of  1500  units. 
The  objective  function  0  is  defined  as 

0  =  ^  cost  coefficient  *  Flow  Rate 
streams  I,  IV,  VIII 

13 


+  35    I        (A.  +  10)°-6  -  ID0"6 
1=1(^7, 8)1 
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TABLE  5 
RESULTS   FuR  THE  PROBLEM  IN   FIGURE  9 


Exchanger  Area 

1  131.50 

2  30.01 

3  413.24 

4  103.24 
2  '  41 . 80 
3*  462.31 
4*  105.06 
l"  131.50 

2  36.10 

3  0.00 
4"  64.67 
5"  0.00 


*'„     =  170,341  F      ,   =  180,878  F         -   50,148 

1  "l  "l 

1  11  n 

T  =  T     =  T     =  190 


0  =    $23451. 25/yr. 
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T  -  350 
F   -   80,000 


Flue   Gas        I 
F-LB-1000 
I   -   700 


ST   -   STKAM,    CU   -    COOL  It  JO    WATER 
LB    -   LCi.'R    BOUND 
I)    -   till).    ALLOWABLE   Al'l'i.OACH   TCMP. 

iat  -  I-'.t:  .".:.■.!. ly  iii-i.a  :a:;.  AiTRuACi  ti:mp. 
*  -  Active  constraints  at  oi-rrn': 


Fig.  10.  An  Example  Exchanger  Network 
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where  A   is  the  area  associated  with  exchanger  i.   In  the  calculation 
of  cost  associated  with  an  exchanger,  the  relation 

cost  =  35[(A  +  10)0,6  -  100,6]  (1) 

is  used  in  preference  to 

cost  =  35  A0"6  (2) 

because  whenever  the  area  of  an  exchanger  currently  at  zero  value,  is 
increased,  the  cost  for  the  exchanger  as  calculated  from  (2)  increases 
abnormally  compared  to  the  change  in  cost  for  the  rest  of  the  system. 
The  modification  as  shown  in  relation  (1)  dampens  this  ill-behaved 
effect. 

There  are  8  decision  variables  for  this  problem  and  the  optimum 
results  after  328  iterations  from  an  infeasible  starting  point.   The 
stopping  criterion  is  a  1  x  10   difference  between  the  worst  and  best 
objective  function  values  in  the  current  set  of  points  retained  by 
the  complex  algorithm.   The  value  of  0  at  the  optimum  is  $152,109/yr. 

Results  for  10  examples  are  shown  in  Table  7.   In  all  the 
examples  the  feasible  point  results  in  very  few  iterations.   It  may  be 
observed  that  the  time  required  for  the  rewriting  of  solution  proce- 
dures after  finding  a  feasible  point  is  relatively  small  as  compared 
to  the  time  taken  for  function  evaluations  during  optimization.   The 
maximum  ratio  of  these  two  times  occur  in  Example  6,  but  it  is  still 
less  than  1/3.   This  observation  indicates  that  the  penalty  paid  for 
rewriting  solution  procedures  whenever  constraint  violations  occur  is 
indeed  very  small. 
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The  figures  shown  for  time  in  Table  7  are  those  required  on 
IBM  360/67.   The  cost  per  second  of  CPU  time  is  about  1.4  cents  and 
the  longest  run  (Example  10)  costs  $8.56  for  complete  execution  while 
Example  1  costs  $0.40.   The  size  limitations  are  50  process  stream, 
25  nodes  and  150  stream  segments  in  the  current  version  of  the  program. 
The  program  is  fairly  well  tested  and  gave  satisfactory  results  when 
used  for  sixteen  different  problems  set  up  by  students  in  a  recent 
design  course. 


CHAPTER  III 
LOCAL  AND  GLOBAL  OPTIMA 

III  .1 .   Statement  of  tlie_  Problem 

In  the  course  of  optimizing  a  chemical  engineering  system 
using  a  conventional  minimum  seeking  algorithm,  one  should  ask  and 
then  attempt  to  discover  if  the  solution  found  is  indeed  a  global 
one.   If  one  is  willing  to  try,  then  a  common  strategy  is  to  re- 
start the  optimization  at  several  different  initial  points,  and, 
if  all  or  most  lead  to  a  single  best  point,  that  point  is  conjectured 
to  be  the  global  optimum. 

A  second  common  strategy  is  to  use  one's  intuition  to  claim 
that  only  one  optimum  is  likely.   This  strategy  is  particularly 
dangerous  when  the  optimum  is  at  the  boundary  of  the  search  region 
where  portions  of  the  system  effectively  disappear  because  the  flow 
through  them  is  zero.   Often  one  models  the  capital  cost  of  a  process 
unit  by  an  equation  of  the  form 

cost  =  C   (Throughput) 

This  form  is  not  convex  and  is  particularly  troublesome  at  zero 

throughput,  where  the  slope  is  infinite  for  cost  versus  throughput. 

If  a  unit  becomes  zero  in  size,  then  a  small  positive  perturbation 

in  its  throughput  has  an  apparent  positive  infinite  effect  on  cost, 

an  effect  usually  large  enough  to  trap  the  optimization  algorithm. 

One  can  of  course  (and  should)  modify  the  cost  equation  to  reduce 
this  problem. 
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These  strategies  may  help  but  still  do  not  guarantee  that  one 
has  found  the  global  optimum.   The  purpose  of  the  work  leading  to  this 
paper  was  to  produce  a  reliable  method  to  determine  if  an  optimal  solu- 
tion is  either  a  global  or  a  local  optimal  solution.   It  was  hoped 
that  the  technique  would  be  computationally  effective,  in  particular 
for  heat  exchanger  networks,  a  class  of  problems  which  commonly  dis- 
plays local  optima. 

The  optimization  problem  relating  to  a  heat-exchanger  network 
is  to  minimize  the  annualized  cost,  0,  of  heat-exchangers  and  utilities 
subject  to  the  approach  temperature  inequality  constraints  and  the 
heat  and  material  balances  representing  equality  constraints. 

The  simple  network  in  Figure  11  (due  to  Grossman  and  Sargent 
(1977))  has  only  one  degree  of  freedom  and  the  network  cost  can  be 
explored  by  choosing  different  values  for  the  temperature  T...   Figure 
12  represents  the  plot  of  costs,  0,  vs.  T„.   At  both  the  points  A  and 
B,  the  Kuhn-Tucker  optimality  conditions  are  satisfied  as  any  slight 
perturbation  away  from  each  of  these  points  results  in  an  increase 
in  0.   However,  the  point  B  clearly  represents  a  local  optimum.   It 
is  quite  likely  therefore  that  an  optimization  algorithm  will  not  find 
the  global  optimum. 

III. 2.   The  Lower  Bound 

To  permit  decomposition  of  a  system  structure,  the  primal  or 

overall  system  optimization  problem  may  be  written  as 

n 
Minimize  0  =  f(q)  =  \    f  .  (q  . ) 
j=l  J   ] 

s.t.  g.  =  x.  -  t.(q.)  =  0   ieO,..   j=l,2,...n 

o1      ±  XX  (j) 

q.ES.    j=l,2,...n         (PI) 
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F   -    10,000 
T   -  600° 


F    -   10,000 
T   -   900° 


F    -    10,000 


<^ — 6- 


C     -   1  I)      -  V.   -   200 

P  12 

Area  -  Q/u  jlT  ,  Cost   ■  35  x  Area0'6 

Cost,    refers    to  exchanger   1. 
Aim:      To  minimize   the   objective    function  i,   where   i  »   Cost     +  Cos 


Fig.   11.    Example  Problem  1 
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Fig,  12.  The  Behavior  of  the  Objective  Function  for 
Example  Problem  1 
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(It  is  important  to  note  that  the  constraints  defining  S   must  be 
complete  enough  to  guarantee  the  boundedness  of  each  subproblem  j.) 

This  representation  of  a  system,  such  as  the  one  in  Figure  13, 
is  based  on  Lasdon  (1970)  ,  but  a  procedure  for  equality  rather  than 
inequality  constraints  is  stressed  as  shown  in  HcGalliard  and 
Westerberg  (1972).   This  stress  follows  from  an  interest  in  decomposed 
system  structures,  as  against,  say,  allocation  of  resources. 

The  problem  is  assumed  to  have  an  optimal  feasible  solution. 
The  Lagrange  function  for  problem  (PI)  is 


L  =  0  -  I  j,         A.[x.  -  t.(q.)] 


j=l  ie0(j) 

which  may  be  rewritten  as 

n  n 

L=      J    ff.(q.)    +        /V        *4t4(q,)    ~         I        AixJ    =      LfMi'h) 


I  [f.CqJ  +     I     \h(^  ~     I     Aixi]  =  .V 

j-l     1     J         ie0(j)  lel(j) 


J     J 


This  Lagrange  function  is  separable  in  q,  and  a  Lagrange  problem  may 
be  defined  which  is  equivalent  to  the  following  set  of  subproblems, 
one  for  each  subsystem  in  the  original  system: 

Minimize   f.(q.,A)     j=l,2,...n 

q.ES.      '   ]  (P2> 

.]   J 

The  solution  to  the  Lagrange  problem  equals  the  sum  of  the  solutions 

of  (P2),  i.e., 

j=l J 
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Fig.    13.     Staged  processes 
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A  dual  function  can  now  be  defined  as 

n   . 
h(A)  =  I    f  (A) 
j  =  l  J 

The  geometric  significance  of  the  dual  h(A)  is  illustrated  in 
Figure  14  for  a  simple  system  like  the  one  shown  in  Figure  13.   The 
dual  h(A)  is  the  intercept  of  the  line  with  slope  A.   Note  that  this 
line  is  a  supporting  hyperplane  at  the  point  (0°,g°).   The  geometric 
significance  of  the  dual  has  been  treated  in  Lasdon  (1970)..  The 
value  of  the  dual  is  always  less  than  or  equal  to  the  value  of  the 
primal  optimum.   Hence  h(A)  can  be  considered  a  lower  bound  on  the 
global  optimum. 

III. 3.   The  Upper  Bound  on  the  Lower  (Dual)  Bound 

Theorem  1 :   If  in  the  space  of  0  vs .  g_,  where  £  is  n- 
dimensional,  a  polytope  (Geoffrion  (1969))  is  formed  in  the  hyper- 
plane passing  through  n+1  support  points,  the  value  of  0  at  the 
intersection  of  this  polytope  with  the  axis  £  =  0  provides  an  upper 
bound  on  the  lower  bound  if  the  point  of  intersection  is  contained 
inside  or  at  the  boundaries  of  the  polytope. 

Discussion:   Figure  15  illustrates  the  ideas  underlying  this 

1       2 
theorem.   g   and  g   are  support  points  to  the  graph  of  0  vs.  g.   The 

1     2 
line  connecting  g   to  g^  is  a  polytope  formed  in  the  hyperplane  pass- 
ing through  these  n+1  =  2  support  points.   This  polytope  or  line 

12 
between  g   and  g   intersects  the  vertical  axis  g  =  0,  and  thus  the 
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PrLnAl    Optical    Valui 


Supporting   hyporplane 
u-ltli   slope   X 


Support   Point 


g  -  vVV 


Fig.    14.    Geometric   Significance   of   the   Dual 
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Upper  Bound  (UB) 


Fig.  15.  Geometric  Significance  of  the  Upper  Bound 
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theorem  says  this  point  of  intersection  represents  an  upper  bound  on 

the  lower  bound. 

1       2 
Proof :   All  the  support  points  of  a  graph  (such  as  g   and  g  ) 

are  contained  on  the  surface  of  a  convex  hull  formed  by  the  support 

hyperplanes  for  this  graph.   Assuming  that  the  polytope  intersects 

the  axis  g  =  0  in  its  interior  or  at  its  boundaries,  a  hyperplane 

intersecting  the  axis  g  =  0  above  this  point  must  intersect  a  part 

of  the  graph,  that  is,  it  cannot  be  a  support  hyperplane.   Thus  every 

support  hyperplane  must  intersect  the  axis  g  =  0  on  or  below  the  point 

the  polytope  intersects  this  axis.   Consequently  the  dual  cannot  exceed 

this  value  of  0  (=VB)    and  UB  is  an   upper  bound  to  the  dual. 

Ill .4 .   Finding  th e  Best  Upper  Bound 

Given  support  points  (0  ,g  )  i=l, 2 , . . . p, pin+1,  the  problem  is 
to  determine  the  following. 

(a)  Find,  if  possible,  a  set  of  n+1  points  that 
yields  an  upper  bound.   Posed  differently,  the 
problem  is  to  find  n+1  points  sue!)  that  the 
polytope  formed  by  them  intersects  the  axis 

g  =  0. 

(b)  If  several  sets  of  n+1  points  qualify  to  pro- 
vide an  upper  bound,  find  the  set  that  yields 
the  lowest  value  for  the  upper  bound. 

Any  point  inside  or  at  the  boundaries  of  the  convex  polytope, 

specifically  the  point  g  =  0,  can  be  obtained  by  a  convex  combination 

of  (n+1)  points  forming  the  polytope.   This  result  has  been  stated  and 

applied  in  Director  et  al.  (1978).   Thus  problem  (a)  can  be  formulated 

as  a  feasible  point  problem  in  linear  programming.   If  a  feasible 

solution  exists  (no  artificial  variables  present  in  the  solution  at 

a  nonzero  level),  then  indeed  the  set  of  n+1  points  does  provide  an 

upper  bound . 
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In  order  Co  solve  problem  (b) ,  if  a  "price"  could  be  associated 
with  each  point  such  that  the  objective  function  reflects  the  value 
of  the  upper  bound  corresponding  to  a  set  of  n+1  points,  then  it  too 
is  a  linear  programming  problem.   The  solution  to  it,  if  one  exists, 
would  yield  solutions  to  both  problems  (a)  and  (b)  simultaneously. 
This  "pricing"  is  indeed  possible. 

For  every  point  on  the  polytope  formed  by  n+1  points,  0  can 
be  represented  by  the  linear  relationship 


T 


0  =  £  a  +  b 


(1) 


where  _a  and  b  are  constants.   The  value  of  0  at  the  intersection  with 
the  axis  ^  =  0  is 


=  b 


(2) 


then 


UB  =  b 


(3) 


assuming  that  g  =  0  lies  inside  or  at  the  boundaries. 

For  n+1  points,  equation  (1)  may  be  rewritten  compactly  as 


0  = 


,n+l 


n+1 


a  +  b 


(4) 


t  _   1'    r  1      n+1,     ,        r   i    .,12      n+1,      _ 
Let  a  =  [a  , . . .a   J  such  that  [a     =  1,  [&  ,£  ,...£   J  a  =  0, 

1=1 


and  a  >  for  all  i.   Premultiplying  equation  (4)  with  a   gives 
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aX0  =  b  (5) 

and  from  equation  (3) 

UB  =  aT0  (6) 

Thus  the  "price"  associated  with  each  point  (01,g1)  is  01. 

The  linear  programming  formulation  to  solve  both  problems  (a) 
and  (b)  simultaneously  is  as  follows: 

Find  u  =  [a  , . . .a  ]  to 


Max  Z  =  -  [0  ,...0F] 


s.t.  [gJ 


=  0 


1     2  r 

a  +  a  +  .  .  .  +  a1 


and  a  ,...,a  ,  >  =  0.   The  upper  bound  UB  =  -  Z. 

The  solution  will  contain  (n+1)  vectors  in  the  basis.   If  a 
feasible  solution  does  not  exist,  artificial  variables  will  be  present 
in  the  solution  at  a  positive  level  and  no  upper  bound  will  exist 
because  the  set  of  points  in  hand  so  far  (g  ,g  , . . .g1)  do  not  surround 
the  origin. 


III.  5.   Improvement  of  the  Upper  Bound 

Let  N  be  the  set  of  the  current  n+1  points  forming  the  poly- 
tope  P,  H  the  hyperplane  passing  through  these  points  with  the  slope 
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_A,  and  UB  the  value  of  the  upper  bound.   The  aim  now  is  to  find  an 
improvement  on  UB . 

If  a  new  point  is  introduced,  it  follows  from  the  simplex 
method  that,  since  the  problem  is  bounded,  a  feasible  solution  and 
hence  a  new  upper  bound  can  be  obtained  by  replacing  one  of  the  points 
in  N  by  the  new  point.   The  question  now  is  how  can  the  new  point  be 
found  so  that  the  new  upper  bound  (UB ' )  is  an  improvement  on  UB, 
i.e. ,  UB'  <  UB. 

Once  again,  the  geometry  of  the  problem  and  the  theory  of 
linear  programming  can  be  used  advantageously.   Assume  for  the  moment 
that  the  solution  to  the  linear  programming  problem  obtained  from  the 
points  in  N  is  not  degenerate.   If  a  support  point  is  determined  for 
slope  A_>  this  new  point  must  lie  on  or  below  H.   This  result  is  ob- 
vious as  a  hyperplane  with  slope  A_  cannot  be  a  support  hyperplane 

above  H.   If  the  coordinates  of  the  new  point  are  (0    ,g    ),  and 

*-u     t     r  n  It         new  .  JA 

the  value  of  0  on  H  at  g_  =   _g_    is  0  ,  then 

0S  >  0n8W 

rtH    rtnew  ,  _  , 

0=0     implies  that  no  further  improvement  can  be  ob- 
tained on  the  upper  bound.   Thus  if  0   is  greater  than  0    ,  it 
follows  that  UB'  will  be  less  than  UB  if  the  new  point  replaces  one 
of  the  points  in  N  by  the  simplex  method.   This  result  can  be 
demonstrated  as  follows: 

Let  the  points  in  N  be  (0  ,g_  )  ,  .  .  .  (0    ,£   ).   The  constraints 


/I   2     n+1.      _  ,.,. 

(£  ,£ £   )  a  =   £  (7) 
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1,2        n+1 

a  +  a.   .  .  .  +  a    =1 


(8) 


Combining  relations  (7)  and  (8) 


G  a 


(9) 


The  new  point  can  be  represented  as  a  linear  combination  of  the  points 
in  the  basis,  that  is, 


„,     n  st   n        .  ~new 

ihe  n+1   element  of  g    gives 


a  ynew 


n+1 
,     r    new 

i  =  l  y  • 
i=i  x 


(10) 


(ii) 


Points  on  H  may  be  expressed  by  equation  (1) 

+  b 


[gT,l] 


0 


For  the  (n+1)  points  in  the  basis,  it  follows  from  equation  (4) 


that 


=  G' 


+  b 


Premultiplying  with  y 


T         T  _, 
new   ,,    new  ^1 
y     0  =  y     G 


Using  equations  (10)  and  (11) 


new   _,      new 
y     0  =  g     ,1 


+  b  I    y\ 
i  =  l  " 


+  b 


(12) 


(13) 


b2 


However,  the  right  hand  side  is  a  point  on  H  at  g     and  thus 

T 
new   M    mH  ,,  ,  N 

y      0  =  0  (14) 


In  the  simplex  method,  the  criterion  for  an  improvement  in  the 

objective  function  by  replacing  a  point  in  the  basis  by  a  new  one  is 

T 
y     0-0    =0-0    >0   (see  Hadley  (1963))  . 

bmce   si      -   v  is   greater    than   zero 

UB'    <   UB 

If  the  linear  programming  solution  corresponding  to  the  points 
in  N  is  degenerate,  that:  is,  the  point  g  =  0  lies  on  one  of  the  bounda- 
ries of  P,  a  difficulty  arises.   In  this  case,  there  is  no  unique 
hyperplane  H  and  consequently  there  may  be  some  degrees  of  freedom 
available  in  calculating  A.   There  is  then  no  guarantee  that  the  new 
point  will  bring  about  an  improvement  in  the  upper  bound.   In  such  a 
case,  several  new  points  may  have  to  be  evaluated  before  an  upper 
bound  with  a  lower  value  is  obtained.   Degenerate  linear  programming 
problem  solutions  occur  in  Example  Problem  3. 

III. 6.   Procedure  to  Investigate  Global  Qptimality 

The  basic  idea  of  this  approach  is  that  the  value  for  the 
global  optimum  must  lie  between  or  at  the  boundary  of  the  upper  and 
lower  bounds.   The  interval  between  the  bounds  is  generally  decreased 
at  every  iteration  by  the  improvement  on  the  upper  bound  and  possibly 
an  improvement  on  the  lower  bound.   If  at  some  point  during  this 
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procedure,  an  optimum  is  found  to  possess  a  value  significantly 
greater  than  the  upper  bound,  an  inference  can  be  made  that  the 
optimum  is  local  and  the  procedure  terminated.   On  the  other  hand, 
if  a  value  for  the  optimum  is  very  close  to  the  lower  bound,  it  may 
be  inferred  that  the  optimum  is  global  and  the  process  terminated. 
A  problem  arises  when  there  is  a  dual  gap  in  the  optimiza- 
tion problem.   It  is  possible  that  the  global  optimum  may  lie  above 
the  upper  bound  at  some  iteration.   However,  the  optimum  is  not 
rejected  as  a  local  optimum  if  it  lies  within  a  certain  interval 
above  the  upper  bound,  say  2%  of  the  value  of  the  optimum.   Previous 
experience  indicates  that  this  modification  is  adequate  for  the  type 
of  problems  considered  in  this  study.   This  behavior  can  be  observed 
in  Example  Problem  2. 

III. 7.   Algorithm 


To  investigate  a  primal  optimal  solution  with  an  objective 
function  0  . 

Step  1    Guess  n+1  different  A's  (A   where  j=l  to  n+1).   Attempt 
to  select  these  AJ  to  obtain  points  g  which  surround  the 
origin. 

Step  2    .For  each  X    ,    evaluate  the  dual_h(AJ),  the  support  point 
g  ,  and  the  objective  function  0J  .   Let  p  =  the  total 
number  of  support  points  (n+1)  in  this  case) .   Find  LB 
where  LB  =  Max  (h(AJ),  j  =  l  to  p)  . 

Step  3     If  0   <  LB  +  e  (c  is  a  small  positive  quantity,  say  .02 
0"),  stop.   0   is  the  global  optimum.   Else  go  to  Step  4. 

Step  4     Formulate  an  L.P.  problem  with  all  the  available  support 
points  and  solve  for  the  optimal  objective  function  (-UB) 
and  the  corresponding  basis  vector.   If  a  feasible  solu- 
tion does  not  exist,  go  to  Step  7.   If  0   >  UB  +  £,  stop. 
The  given  solution  is  a  local  optimum.   If  the  L.P.  solu- 
tion is  degenerate,  go  to  Step  6. 
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Step  5     In  the  L.P.  solution,  let  the  vectors  in  the  basis  be 

g  ,g  ,g  , . . . g    .   Set  up  n  simultaneous  equations  in 
unknowns  Ap+I  as  follows 

,    2    l.T  .p+1  rt2    A 

(g   -  g  )   A     =  0-0 


(15) 


,3    l.T  ,p+l    ,A3    A 
(g   -  g  )   X1    =0-0 

,  4    l.T  .p+1    ,.4    -1 
(g   -  g  )  Y        =0-0 

,  n+1     l.T  ,p+l    „,n+l    A. 
(g    -  g  )   AF   =0    -  0 

Solve  for  A'   .   Find  the  dual  h(A    ),  the  support 
point  gP   ,  and  the  objective  function  0P   .   Set 
p  =  p  +  1.   Find  LB  where 

LB  =  Max  (h(A_J)  ,   j=l,2,  .  .  .p) 

Go  to  Step  3. 

Step  6    Select  all  vectors  from  the  basis  which  are  present 

in  the  L.P.  solution  in  Step  4  at  a  nonzero  level.   Let 
these  be  points  g',g~,...gm.   Since  only  m  (<  n+1)  such 
points  exist,  the  AP+J-  to  be  chosen  can  come  anywhere 
from  within  an  n+1  -  m  dimensional  subspace.   Either 
initiate  (if  first  time  through  this  step  with  above 
nonzero  basis  vectors)  or  continue  a  systematic  search 
over  the  subspace  to  find  the  next  A's  to  use.   A  A 
from  the  appropriate  subspace  is  found  by  having  it 
satisfy 

,2    l.T  ,  p+1    rt2    .1 
(g   -  g  )   Ap    =0-0 

(g3  -  gV  AP+1  -  03  -  01  (16) 

(gm  -  gV  AP+L  =  01"  -  01 

plus  any  additional  n  +  1  -  m  independent  specifications. 
Go  to  Step  3. 

Step  7     Select  a  set  of  A     to  find  a  new  support  point  g    , 
the  dual  h(AP+1),  and  the  objective  function  0P+1 .   Set 
p  =  p  +  1.   Let  LB  =  Max  (h(A.I),  where  j=l  to  p)  .   Go  to 
Step  3.   (Note  this  step  is  actually  the  same  as  Step  6 
but  with  no  constraints  of  form  (16)  being  written,  i.e., 
one  should  search  systematically  over  the  entire  space 
for  the  next  A's.) 
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I I I. 8.   Examples 

Example  Problem  1 

Figure  11  represents  the  Example  Problem  1,  and  its  subsystems 
are  shown  in  Figure  16. 

h(A)  =  Min  (cost,  +  Ay)  +  Min  (cost,  -  Ax) 


where  A   is  the  price  associated  with  the  variable  T  . 

81  =  (x2  -  yl} 
On  optimization  of  the  original  system  the  two  possible 
solutions  are 


Case:   1   T.  =  400°    T2  =  900°    T3  =  300°    0   =  286.93 
Case:   2   T  =  600°    1'2  =  700°    T  =  100°    0*  =  189.3 

Apply  the  algorithm  to  Case  1,0   =  286.93 
Step  1   A   =  -4.14  ;    A   =  4.14   and   p  =  2  points. 
Step  2   gj  =  -200      h(xh  =  -224.6   01  =  476.2 

g^  =  +200      h(A^)  =  -882.0  02  =  0 

LB  =  -224.6 

Step  3   0*  >  LB  +  e.   (Let  c  =  0.02  0*  =  5.74.) 

Step  4   UB  =  238.1 
0  >  UB  +  e 


Hence  the  solution  in  Case  1  is  a  local  optimum. 


6b 


F    »   10,000 

r  -  6ou° 


F    -    10,000 

r  -  300° 


<> 


F   -    10,000 


~ 0" 


T, 


100  s    y     S    300 


100  s    x     s    J10 


FiR.  1*j.  Subsystems  for  Example  Problem  1 
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Apply  the  algorithm  to  Case  2,  0     189.3  Steps  1,  2  and  3  are 
the  same  as  in  Case  1. 

Step  4   UB  =  238.  1 
0   <  UB  +  E 

Step  5   A3  =  -1.19 

h(X?)  =189.3   ,    gj  =  0   ,    03  =  189.3 
P  =  3 
LB  =  189.3 
Step  3  0  =  189.3  <  LB  +  e 

Hence  the  solution  in  Case  2  is  a  global  optimum. 

Example  Problem  2 

Figure  17  represents  this  problem  and  the  subsystems  are  shown 
in  Figure  18. 

h(A)  =  Min  (Area   +  Ay)  +  Min   (Area„  +  A9y    -  A  x  ;) 
yll  yi2,X22 

4-  Min  (Area.,  -  Ax,,,,) 
X23 


A   is  the  price  associated  with  T   and  A   is  the  price  associated  wi 


th 


T  . 
2 


81  =  (x22  ~  yll)   and   82  =  (x23  "  y12^ 


On  optimization 


Tx   =    181.9°   T2  =  295°    t   =  218.1°   t?  =  286.4° 
t3  =  395.5°   and   0   =  7049. 
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T  =   100"  /^-\ 


-u — - i 


tcJ 


FC     -   10      tor  all   the   stream 
P 

U,    -   120  U„   -  80  IL    ■  40 


Arua  refers  to  exchanger  1 
Area  -  Q/D  il^  ,  d  -  Are^ 
Aira:         To  minimise   4 


>.rc»,  +  Area, 


[•'ig.    17.     Exsrunle  Problem 
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Subsystem  1 

T  >  300° 


Subsystem  2 

T  -  «iO° 


Subsystem  3 

I  -  600° 


^A-       *_A 


100  £  y    s  300 
100°  St,   s  300° 


100  £  x„„  *  300 
22 

100  £  y   £  400° 

100°  s  t,  £  400° 


100  £  x   £  400 
200°  st,  £  600° 


Fig.  18.  Subsystems  for  Example  Problem  2 
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Application  of  the  Algorithm:   0   =  7049 


Step  1 

A1  = 

-10 
0 

2 

0 
10 

A3  = 

-15 
_-30_ 

Step  2 

1 

S   = 

200 
100 

2 
S   = 

0 
300 

3 
1  = 

-94.6 
-188 

h(A  )  =  500 
01  =  2500 


Step  3   0   >  LB  +  e 
Step  4  UB  =  9001.46 

■k 

0  <  UB  +  e 


h(A2)  =  -500  h(A3)  =  5787 


Step  5   A4  = 


-21.67 
-21.67 


0  =  2500 
LB  =  5787 
Let  z   =   140 


■112.3 
132 


=  12,846 


h(AH)  =  5585   0  =  5158     p  =  4 
LB  =  5787 
Step  3   0   >  LB  +  e 
Step  4   UB  =  7173.34 

0   <  UB  +  e 

13      4 
Step  5   The  points  in  the  basis  are  g  ,  g   and  g  . 


-11.04 
-24.65 


122.87 
7J 


h(A5)  =  6640.3   05=  3533.67   p  =  5 


LB  =  6640.3 


P  =  3 


Step  3   0   >  LB  +  E 


71 


Step  4   UB  =  6928.94 
0   <  UB  +  e 


3   4      5 
Step  5   The  points  in  the  basis  are  g  ,  g   and  g  . 


-13.34 
-24.76 


111.79 
71 


h(A6)  =  6916.26   06  =  3668.11 


P  =  6 


LB  =  6916.26 

Step  3   0   <  LB  +  e 

Hence  the  solution  represents  a  global  optimum. 

Example  Problem  3 

Figure  2  represents  this  problem  and  the  subsystems  are 
shown  in  Figure  19. 

y   ,  x   and  price  A   are  associated  with  F? 

y,n,  x„,  and  price  A„  are  associated  with  (F.  x  T , ) 
712   24  2  4    4 

y„9,  x„s  and  price  A~  are  associated  with  T  ? 

y,.-,,  x,.,  and  price  A,  are  associated  with  T,  c 
1326  4  1j 


h(A)  =  Min   (Cost   -  ^ix?o  +  ^o^l  2  +  ^-^22^  + 

X22'y12'y22 

Min   (Cost„  +  Ay  -.  -  A  x„,  +  A,  y  „)  + 


'll>A24"l3 

Min   (Cost-  +  0.4F    -  ^X^)  + 


~25 

Min   (Cost.  +  0.2FTTT  -  A.x„,)  + 
4        III     4  2b 


^26 


and 
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?    •   lO.OL'O 

r  -  200° 


Subsystem   2 


yil 

T 

X24 

F  -  10,000 

1 

10,000 

200° 

* 

T  -  800 

T 

■  40U^ 

2500  <-    *22'yi2   s    "00 


300      £   y        S    500° 


2500  i   y1,,x  4  s    7500 
300°  s   y       s    500° 


Sub  ay  stem  t* 

F   -   10,000 


300     s    x„c   s    500 


Fig.  19.  Subsystems  for  Example  Problem  3 
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gl  =  (x22  ~  yll} 

82  =  (x24  "  y12} 

g3  =  (x25  -  y22) 

84  =  (x26  "  y13} 


On  optimization  there  are  two  possible  solutions: 


Case:    1   F„  =  2500   F0  =  7500   T,  =  Tc  =  400° 
2  J  4     5 


12 


=  300° 


T15  =  500°   Fn  =  2500   Fm 


4154 


Case ; 


2   Y     =    7350   F3 


2650   T       =  500' 


T15  =  300c 


F1I  =  ° 


FIII  =  25°°  T5  =  /(23° 


T,   =  392° 

4 


3538 


Application  of  the  Algorithm  to  Case  1    0     4154 


Step    1      A      = 


0 

0 

-10 

•10 


A2   = 


0 
0.000555 
0 
0 


0 

0 

10 

0 


A4   = 


-0.1 

0.001 

-10 

-10 

A5   = 


0 
0.001 
-10 
-    5 


P   =    5 


Step    2      gJ 


5000 

-2,000,000 
■200 
■200 


2 
g      = 


-5000 

2,000,000 

200 

200 


2345 
-1,876,000 
■200 

200 


■5000 
2,000,000 

0 

0 


-2345 

1,876,000 
■200 

200 
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hCA1)  =  2350          01  =  6350 

h(A2)  =  230           02  =  1340 

h(A3)  =  1917          03  =  3917 

LB  =  2350 
h(A4)  =  2118         0  =  4618 

h(A5)  =  2188          05  =  5064 

Step  3 

0"  >  LB  +  e     Let  E  =  80 

Step  4 

UB  =  3845 

0   >  UB  +  e 

Hence  Case  1  corresponds  to  a  local  optimum. 

* 

Application  of  the  Algorithm  to  Case  2   0    3538 

Steps  1  to  3  same  as  in  Case  1 . 

Step  4 

UB  =  3845 

0   <  UB  +  e 

LP  solution  is  degenerate 

Step  6 

Nonzero  points  in  the  basis  of  LP  solution  are  g   and  g  . 

Let  the  n+1  points  for  finding  A   be  g  ,  g  ,  g  ,  g  ,  and  g  . 

A6  = 

-0.065 

6 
S   = 

0 

0.000225 

0 

-9.81 
-6.58 

-200 
200 

h(A6)  =  3338           0  =  4037 

p  =  6                  LB  =  3338 

Step  3 

0   >  LB  +  e 

Step  4 

UB  =  : 

5845 

0  <  ub  +  e 

LP  solution  is  degenerate 
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1      2 
Step  6   Nonzero  points  in  the  basis  of  the  LP  solution  are  g   and  g 

Let  the  n+1  points  for  finding  A   be  g  ,  g  ,  g  ,  g  ,  and  g  . 


A7  = 

0.128 

0.000709 

-8.68 

-7.72 

7 
8  = 


-2345 

4690 

0 

0 


h(A7)  =  3534 


P  =  7 


0  =  4564 
LB  =  3534 
Step  3   0   <  LB  +  e. 

Hence  Case  2  represents  the  global  optimum. 

III. 9.   Discussion 

In  the  first  and  the  third  examples,  the  local  optima  are 
recognized  very  quickly.   The  global  optimum  in  each  of  these 
examples  is  confirmed  albeit  with  a  greater  number  of  steps.   In  the 
second  example,  a  dual  gap  (global  optimum  =  7049,  upper  bound  =  6928.9) 
is  present.   However,  the  algorithm  proves  effective.   It  may  also  be 
noticed  in  the  second  example  that  an  improvement  results  in  the  upper 
bound  at  every  iteration.   This  is  guaranteed  in  the  absence  of  de- 
generacy.  Degenerate  solutions  occur  in  the  third  example  and  several 
points  may  have  to  be  evaluated  before  obtaining  an  improvement  in  the 
upper  bound. 


CHAPTER  IV 

PROCESS  SYNTHESIS  USING  STRUCTURAL  PARAMETERS: 

A  PROBLEM  WITH  INEQUALITY  CONSTRAINTS 

In  this  chapter  a  problem  with  the  structural  parameter  method 
for  process  synthesis  shall  be  discussed.   The  problem  is  almost 
obvious  once  stated,  but  it  has  not  been  emphasized  in  the  literature. 
It  can  seriously  affect  the  expected  results.   The  discussion  here  does 
not  imply  that  the  method  is  valueless,  it  simply  stresses  that  care 
must  be  taken  to  ensure  that  the  method  is  really  useful  for  a  given 
problem. 

A  system  of  interconnected  process  subsystems  can  be  modeled 
using  structural  parameters  which  are  defined  by  the  following 

equations 

N 
x.  =   y   a. .x!  i=l,2,3, . . .N 

J  =  l 


£  a. .  £  1  ,    y  a. .  =  1 

1=1 


(1) 


where  x.  and  x!  are,  respectively,  the  input  and  output  variables  of 
the  i-th  subsystem,  N  is  the  total  number  of  subsystems  in  the  entire 
system  and  the  parameters  a.,  are  the  structural  parameters;  that  is, 
each  is  the  fraction  of  the  output  stream  of  the  j-th  subsystem  which 
flows  into  the  i-th  subsystem. 

By  means  of  structural  parameters,  a  system  synthesis  problem  can 
be  transformed  into  a  non-linear  programming  problem  with  continuous 
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decision  variables.   This  idea  has  been  stated  and  demonstrated  in 

the  studies  by  Umeda,  Hirai  and  Ichikawa  (1972),  lchikawa  and  Fan 

(1973),  Osakada  and  Fan  (1973),  Mishra,  Fan  and  Erickson  (1973),  and 

Himmelblau  (1975).   In  order  to  obtain  an  optimal  structure,  redundant 

subsystems  are  usually  inserted  into  a  "super"  structure  in  which  it 

is  hoped  that  the  optimal  structure  is  imbedded.   For  example,  Figure 

20  illustrates  one  use  of  the  method  for  synthesizing  a  heat  recovery 

network.   The  problem  data  and  the  stream  specifications  are  shown 

in  Table  8,  and  the  problem  is  described  in  the  tradition  of  problems 

like  4SP1  in  Masso  and  Rudd  (1969)  and  Lee  et  al.  (1970).   A  decision 

is  to  be  made  on  whether  or  not  to  use  cold  stream  II  to  aid  in 

cooling  hot  stream  V.   The  problem  is  formulated  with  potentially 

redundant  exchangers  2,  4  and  5  and  the  structural  parameter  a   . 

a0/.  is  the  fraction  of  cold  stream  II  that  goes  through  exchanger 
26 

2.   The  strategy  is,  in  this  formulation,  to  convert  the  discrete 
decision  of  whether  or  not  to  use  cold  stream  II  to  aid  in  cooling 
hot  stream  V  into  a  problem  of  optimizing  over  the  continuous  deci- 
sion variable  a„, .   To  avoid  a  cost  discontinuity  caused  by  intro- 
Zb 

ducing  an  exchanger  we  must,  of  course,  require  the  cost  of  an 
exchanger  to  approach  zero  as  its  area  does. 

The  notion  that  we  have  converted  a  discrete  decision  problem 
to  a  continuous  one  is  invalid  here.   "It  may  be  observed  that  hot 
stream  V  has  sufficient  heat  content  to  drive  cold  stream  I  to  its 
final  temperature.   However,  the  transfer  of  heat  in  exchanger  I  is 
limited  because  an  inequality  constraint  in  exchanger  2  has  to  be 
satisfied.   This  constraint  is  the  approach  temperature  requirement 
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II,  Cold 
t  -  504.00 
T  -  311.11 


Cost   -   16,945   S/yx 


VII,    Heating   Utility 
T  -  588.89 


Fig.   20.    The  Optimal  Values  of  the  Parameters 
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TABLE  8 
STREAM  SPECIFICATIONS  AND  PROBLEM  DATA 


'^J'-rimoN 

.  is; is .... 

STREAM 

flow 

g/» 

I 

II 

III 

IV 

V 

vr 

¥11 

1260. 00 

504.00 

1260.00 

Unknown 

2520.00 

Unknown 

Unknown 

Iclet  Teapeiature 

i 

311.  a 

311.11 

422.22 

311.11 

644.  U. 

644.44 

588.3-7 

Outlet  Traperatire 

i 

588.89 

533.33 

477.78 

355.56 

477.78 

644.44 

588.89 

Bolllnj  Point 

K 

755.56 

755.56 

755.56 

373.33 

755.56 

644.  44 

588.89 

Liquid  Heat  Capacity 

kJAg  s 

A. 19 

4.19 

4.19 

4.19 

4.19 

4.19 

4.19 

Film  lk.it  Tr.w-.3fer 
Coefficient 

W/m2  K 

1703.49 

1703.49 

1703.49 

1703.49 

1703 .49 

8517.45 

8517.45 

Coot 

8 

0.00 

0.00 

o.oo 

l.lOxlo"7 

0.00 

22.05*10~7 

5.20xlO~7 

Heat  of  Vaporization 

kJAg 

6?/.  80 

697.60 

697.80 

1786.368 

697.80 

1628.20 

1395.60 

Inlet  Vapor  fraction 

0.00 

0.00 

0.00 

o.oo 

0.00 

1.00 

1.00 

Dutlet  Vapor  Fraction 
i —  _ 

0.00 

0.00 

0.00 

0.00 

0.00 

o.co 

O.CO 

Approach  Temperature     =     2.78  K 

Beat  exchanger  coot  equation     =        (a  A    ) 

where         =  arum-vl   rata  of  return  =  0.1 

a     =  350 

b     =  0.6 

Equijwiont  down  tins  =  280  hr./yr. 
The  cost  of  the  network  is  the  cost  of  utility  streams  +■   tho  cost  of  the  exchangers,   In  l/yr. 
Ill   heat  archangel's  are  ausuaed   to  bo  counter-current 
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between  the  inlet  hot  stream  and  the  outlet  cold  stream  (T   > 

533.33°  +  D,  where  D  is  a  minimum  allowed  approach  temperature  of,  say, 

2.78K).   By  numerical  computation  one  finds  that  the  overall  cost  can 

be  lowered  if  a        is  increased  (see  Figure  21).   This  increase  in  a„, 
£D  26 

lowers  the  flows  of  cooling  utilities  IV  and  VII.   Ultimately,  when 

the  flow  of  cooling  utility  IV  is  zero,  an  optimal  structure  is  obtained 

with  a  cost  of  $16,945/yr.   At  this  point  a_,  =  0.688. 

26 

Note,  however,  that  when  the  flow  of  the  cold  stream  through 
exchanger  2  is  zero,  and  the  heat  exchanger  2  is  totally  neglected, 
a  network  with  a  cost  of  $6864/yr.  is  obtained  as  shown  in  Figure  22. 
Thus,  a  significant  discontinuity,  which  was  hoped  to  have  been 
eliminated,  has  reappeared.   It  should  be  emphasized  that  the  data 
for  the  problem  and  its  formulation  are  important  only  in  that  they 
are  plausible  and  that  they  help  make  this  point. 

The  main  observation  may  be  summarized  as  follows.   When  certain 
subsystems  are  rendered  redundant  during  optimization  (by  the  asso- 
ciated flow  or  a  split  fraction  taking  on  a  value  of  zero)  and  if 
inequality  constraints  are  associated  with  these  subsystems  which 
force  constrained  behavior  elsewhere,  discontinuities  very  likely 
still  exist  in  the  problem.   One  must  still  make  a  discrete  decision 
about  whether  or  not  to  introduce  that  subsystem.   This  observation 
means  that,  if  the  structural  parameters  are  to  be  used  for  a  syn- 
thesis problem,  and,  if  inequality  constraints  are  involved,  the 
problem  has  to  be  formulated  very  carefully,  if  indeed  it  can  be, 
to  make  it  a  continuous  one. 
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Fig.  21.     The  Discontinuity  in  Optimization 
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Fig.  22.  The  Optimal  Network  When  Ignoring  the  Approach 
Temperature  Requirement  at  Null  Exchanger  2  of 
Figure  20 


CHAPTER  V 
CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FURTHER  RESEARCH 

The  optimization  strategy  chosen  for  EROS  proves  to  be  efficient. 
It  is  the  author's  belief  that  the  number  of  steps  required  for  con- 
vergence is  significantly  lowered  by  rederiving  a  solution  procedure 
every  time  a  constraint  violation  occurs,  and  by  the  use  of  'restriction' 
(Geoff rion  (1970)).   A  very  useful  feature  in  EROS  is  its  ability  to 
find  a  feasible  starting  point,  if  none  such  is  provided  by  the  user. 
The  solution  yielded  by  EROS  can  account  for  portions  of  the  network 
that  already  exist  and  for  irregularities  likely  to  occur  in  the 
process  streams.   Considering  the  nature  of  the  problem  treated,  cost 
per  a  typical  run  of  EROS  seems  small.   The  use  of  program  may  also 
be  extended  in  carrying  out  synthesis  via  structural  parameters  but 
the  problem  must  be  formulated  very  carefully  as  mentioned  in  Chapter 
IV.   In  the  results  shown  in  Table  7,  a  feasible  starting  point  is 
obtained  in  very  few  iterations.   A  feasible  starting  point  may  also 
be  provided  by  the  user  as  in  example  10  of  Table  7.   As  indicated 
in  Chapter  II,  the  time  required  for  the  rewriting  of  solution 
procedures  after  finding  a  feasible  point  is  relatively  small  as 
compared  to  the  function  evaluations  during  optimization.   Hence, 
the  equation  solving  approach  followed  is  fully  justified. 
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The  practical  applicability  of  a  program  such  as  EROS  can  only 
be  considered  complete  if  a  capability  of  handling  power  generation 
and  power  intensive  units  is  incorporated.   Hence  there  is  a  need  for 
adding  unit  models  such  as  turbines  and  compressors.   Hereafter  the 
pressure  changes  will  be  important  and  a  pressure  variable  must  be 
associated  with  each  segment.   These  changes  might  also  call  for  a 
better  modeling  of  a  heat-exchanger  unit  in  which  the  heat-transfer 
coefficients  may  no  longer  be  evaluated  by  simplistic  relations. 
Currently  the  cooling  curves  for  streams  are  approximated  by  three 
linear  regions.   A  more  accurate  description  could  be  provided  for 
streams  by  specifying  enthalpy  versus  temperature  information  for  the 
range  under  consideration.   A  linear  interpolation  may  be  assumed 
between  two  adjacent  points  provided  in  this  specification. 

In  Chapter  III  an  algorithm  is  presented  which  can  effectively 
identify  if  an  optimum  point  is  a  global  optimum  or  a  local  one.   It 
is  practical  as  given  only  when  applied  to  systems  of  interconnected 
subsystems,  fortunately,  a  problem  form  which  is  common  in  engineer- 
ing design.   The  algorithm  develops  a  lower  (dual)  bound  using  the 
ideas  of  Brosilow  and  Lasdon  (1965)  wherein  the  total  system  is  de- 
composed into  subsystems,  each  of  which  must  be  globally  optimized. 
If  the  subsystems  are  small  enough,  this  requirement  is  readily  satis- 
fied.  Frequently  the  optimal  lower  (dual)  bound  is  equal  to  the 
optimal  upper  (primal)  bound,  but  it  also  often  lies  strictly  below 
because  of  nonconvexities  in  the  problem.   Previous  experience  on 
actual  engineering  problems  indicates  that  this  difference,  caused  by 
a  duality  gap,  is  small,  i.e.,  only  a  small  percent  of  the  system  cost. 
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The  algorithm  also  develops  an  upper  bound  to  this  lower  bound. 
Subject  to  the  tolerances  required  because  of  the  duality  gap,  an 
optimum  more  than  a  few  percent  above  this  upper  bound  is  declared 
to  be  a  local  optimum,  and  one  at  or  just  above  the  lower  or  dual 
bound  is  considered  to  be  a  global  optimum.   If  no  discrimination  is 
possible  at  the  current  step,  the  algorithm  provides  a  means  to 
guarantee  improvement  of  the  upper  bound  for  nondegenerate  problems, 
the  lower  bound  usually  being  improved  at  the  same  time.   The  three 
example  problems  demonstrate  that  the  algorithm  is  surprisingly 
effective. 

Relating  to  the  algorithm  presented  in  Chapter  III,  a  study  is 
recommended  to  discover  the  global  optimum  having  discovered  that 
the  given  primal  optimum  is  local.   It  would  be  desirable  to  predict 
a  value  of  A  yielding  the  global  optimum  by  perceiving  a  trend  in 
the  value  of  X    from  the  information  gathered. 
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APPENDIX  A 
DESCRIPTION  OF  THE  INPUT  AND  OUTPUT  FEATURES 


DATA  INPUT  FOR  THE  PROGRAM  EROS 


Card 


Variable   Card  Format   Description 


First  Card 


FL 


SC 

Next  Card   MNODES 
MSEGS 

NS 
LS 

NESTK 


NSTOP 


4F8.1      Scaling  factor  for  Flow  Recommended 
FL  =  Minimum  Flow 

Scaling  factor  for  Temperature 
Recommended  T  =  Minimum  Temperature 

Scaling  factor  for  split  fraction 
Recommended  A  =  0.2 

Scaling  factor  for  slack  variables 
Recommended  SC  =  Minimum  Temperature 

11  14      Number  of  nodes 

Number  of  stream  segments 

Number  of  streams 

Number  of  information  items  in  a 
stream  =  18 

(length  of  stream  vector)  NESTK 
=  0  if  a  starting  point  is  not 
provided 

=  number  of  constraints  specified 
for  a  starting  point 

=  0  if  the  program  is  to  be  executed 
until  completion 

=  1  if  the  program  is  to  be  stopped 
after  one  iteration  -  to  permit 
a  check  that  the  input  has  been 
read  correctly 
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Card 


Variable   Card  Format   Description 


NSEGS 


Number  of  segments  (note:   not 
streams)  on  which  an  upper  and  a 
lower  bound  on  enthalpy  is  to  be 
specified 


INF 


1NH 


Total  number  of  entries  in  the 
vector  KNF  (see  the  description 
of  KNF) . 

=  0  if  the  option  is  not  used 

Total  number  of  entries  in  the 
vector  KNH  (see  the  description 
of  KNH) . 

=  0  if  the  option  is  not  used 


Optional 

Cards 
(when 
NESTK  i 
0) 


1NA 


JWRIT 


KESTK  (I)  918 
1=1,  NESTK 


Total  number  of  entries  in  the 
vector  KNA  (see  the  description 

of  KNA) . 

=  0  if  the  option  is  not  used 

=  0  if  a  complete  output  is  to 
be  printed 

=  1  if  the  output  during  optimi- 
zation is  to  be  suppressed. 

List  of  codes  for  the  initial 
constraints  to  be  held  (see  the 
section  on  how  to  "code"  a  con- 
straint) 


Next  set 
of  Cards 
=  3xNS 
where  NS 
is  the 
number  of 
streams 


STPAR  (I, J)  8F8.1 
1=1,  18 
J=l,  NS 


Stream  Information 
A  set  of  3  cards  is  necessary  for 
every  stream,  with  up  to  8  entries 
per  card 


For  Any  Stream  .J 


1st  Card 


STPAR  (1,J) 


-2.  for  undefined  flow  of  a 
process  stream 

-1 .  for  undefined  flow  of  a 
utility  (steam  or  cooling  water) 
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Card 


Variable   Card  Format   Description 


2nd  Card 


STPAR  (2, J) 

STPAR  (3, J) 

STPAR  (4, J) 

STPAR  (5, J) 

STPAR  (6, J) 

STPAR  (7, J) 

STPAR  (8, J) 

STPAR  (9, J) 

STPAR  (10, J) 

STPAR  (11, J) 

STPAR  (12, J) 

STPAR  (13, J) 

STPAR  (14, J) 


The  undefined  flow  is  identified 
in  two  different  ways  because  the 
costing  of  the  streams  is  different 
in  each  case. 

The  cost  of  a  process  stream  need 
be  found  once  regardless  of  the 
number  of  units  it  passes  through. 
A  utility  stream  cost  has  to  be 
determined  at  every  heater  or 
cooler  that  exists  in  which  it  is 
a  utility. 

Inlet  temperature 
=  -1.  if  undefined 

Output  temperature 
=  -1.  if  undefined 

=  0.   Inlet  pressure 
(Not  used  at  present) . 

=  0.   Outlet  pressure 
(Not  used  at  present). 

Inlet  vapor  fraction 
=  -1.  if  the  inlet  temperature 
is  undefined 

Outlet  vapor  fraction 
=  -1.  if  the  outlet  temperature 
is  undefined 

Liquid  heat  capacity 

Vapor  heat  capacity 

Heat  of  vaporization 

Dew  point  temperature 

(Must  be  at  least  1°  greater  than 

the  bubble  point  temperature) 

Bubble  point  temperature 

Liquid-side  film  heat  transfer 
coefficient 

Vapor-side  film  heat  transfer 
coefficient 
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Card 


Variable     Card  Format   Description 


3rd  Card 


STPAR  (15, J) 
STPAR  (16, J) 
STPAR  (17, J) 


Next 
Card 


STPAR  (18, J) 


Al 


4F16.8 


A  2 


Two  phase  film  heat  transfer 
coefficient 

Cost  coefficient  ($/lb.*  no.  of 
operating  hrs.  per  yr.) 

Lower  bound  on  flow.   The  lower 
bound  must  be  set  to  quantity 
greater  than  zero  (say  1%  of  total 
flow)  if  the  stream  is  split. 
This  helps  avoid  numerical  dif- 
ficulties . 

Upper  bound  on  flows 
(Can  be  very  large) 

Step-size  for  Kuhn-Tucker  test. 
Recommended  value  =  0.01  (change 
in  scaled  values  for  variables) . 

Initial  step-size  in  scaled 
variable  values  for  the  complex 
method.   Recommended  value  =0.1 


A3 


A  4 


Minimum  allowable  approach 
temperature  in  degrees 

Stopping  criterion  for  the 
optimization  (change  in  objec- 
tive function) .   Recommended 
value  =  0.00001 


Next 
Card 


CA 


3F16.8  Cost  multiplier  for  an  exchanger 
(heater  or  cooler)  x  annual  rate 
of  return. 


EXPO 


Exponent  for  the  calculation  of 
the  cost  of  an  exchanger,  heater 
or  cooler.  Recommended  value  = 
0.6. 


CB 


Modification  in  the  cost  function 
of  an  exchanger,  heater  or  cooler. 
Recommended  value  =  10.   (Aids  in 
avoiding  the  creation  of  local 
minima. ) 

Cost  of  an  exchanger,  heater  or 
cooler  =  CA*((Area  +  CB)**0.6  - 
CB**0.6) 
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Card 


Variable 


Card  Format   Descrii 


Next  Set  of   ARK(I) 
Cards        I  =  1.MN0DES 


Optional 
Cards 

(when  INA  / 
0) 


6F10.2      This  vector  contains  the  infor- 
mation whether  area  is  specified 
for  a  given  node. 

=  0.  if  area  is  not  specified 
=  value  of  area,  if  specified 
6  area  entries  per  card.  For 
mixers  and  splitters,  a  value 
of  zero  is  specified. 


KNA(I)         1114        Gives  information  about  ex- 
I  =  1,INA  changers,  heaters  and  coolers 

that  are  to  be  related  by  areas 

Structure  of  KNA  vector 
Entry  1 

Item  1:   Number  of  nodes  in 
entry    (say,  n,) 

Items  2,3,4,  .  .  .rij  +  1 

Node  numbers  in  ascending  order. 
Repeat  for  entry  2  etc.   New 
entries  do  not  start  a  new  card. 


Optional     KNF(l) 
Cards        I  =  l.INF 
(when  INF  $ 
0) 


Optional     KNH(I) 
Cards        I  =  l.INH 
(when  INH  t 
0) 


1114 


Gives  information  about  un- 
related segments  which  are  to 
have  equal  flows.   No  segment 
flows  are  to  be  listed  here  as 
equal  if  they  are  already  equal 
by  the  nature  of  the  flowsheet. 

Structure  of  KNF  vector 

Entry  1 

Item  1:   Number  of  stream  seg- 
ments in  the  entry  (say,  n-i) 

Item  2,3,  ..  .n  +  1 

Stream  segment  numbers  in 
ascending  order. 

Repeat  for  entry  2  etc.   New 
entries  do  not  start  a  new  card. 

The  description  of  this  vector 
is  the  same  as  that  of  KNF  with 
the  exception  that  enthalpies 
are  involved  instead  of  flows. 
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Card 


Variable      Card  Format   Description 


No.  of  Cards  NARRAY(I,J)     1114 

=  MNODES  J  =  1,11 

(No.  of  1=1,  MNODES 

nodes) 


One  card  for  each  node. 
Every  card  has  11  entries. 


For  node  I  which  represents  an  exchanger  heater  or  cooler 
N ARRAY (1,1) 


NARRAY(I,2) 
NARRAY(I,3) 

NARRAY (I, 4) 
NARRAY(I,5) 

NARRAY(1,6) 
NARRAY (1,7) 

N ARRAY (I, 8) 
NARRAY (I, 9) 

NARRAY (I, 10) 
NARRAY (1,11) 


Stream  I.D.  for  the  hot  stream 
segments 

hot  input  segment  number 

source  node  for  hot  input 
segment 

hot  output  segment  number 

destination  node  for  hot  output 
segment 

cold  input  segment  number 

source  node  for  cold  input 
segment 

cold  output  segment  number 

destination  node  for  cold 
output  segment 

Type  of  node 

=  1  for  a  heat  exchanger 

=  2  for  a  heater 

=  3  for  a  cooler 

Stream  I.D.  for  the  cold  stream 
segments 


For  node  I  which  represents  a  splitter 
NARRAY (1,1) 


NARRAY (I, 2) 


Stream  I.D.  for  the  hot  stream 

segment 

=  0  for  a  cold  stream 

hot  input  segment  number 
=  -1  for  a  cold  stream 
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Card       Variable      Card  Format     Description 

NARRAY(I,3)  source  node  for  hot  input  segment 

=  -1  for  a  cold  stream 

NARRAY(I,4)  the  first  output  segment  number 

NARRAY(I,5)  destination  node  for  the  first 

output  segment 

NARRAY(I,6)  cold  input  segment  number 

=  -1  for  a  hot  stream 

NARRAY(I,7)  source  node  for  cold  input  seg- 

ment 
=  -1  for  a  hot  stream 

NARRAY(I,8)  the  second  output  segment  number 

NARRAY(I,9)  destination  node  for  the  second 

output  segment 

NARRAY(I,10)  Type  of  node 

=  4  for  a  splitter 

NARRAY(I,11)  Stream  I.D.  for  the  cold  stream 

segments 
=  0  for  a  hot  stream 

Note:   the  order  for  the  output 
segments  can  be  arbitrary. 

For  node  I  representing  a  mixer 

NARRAY(I,1)  Stream  I.D.  for  the  hot  stream 

segments 
=  0  for  a  cold  stream 


NARRAY(I,2) 
NARRAY(I,3) 

NARRAY(I,4) 

NARRAY(I,5) 

N ARRAY (I, 6) 
NARRAY(I,7) 


the  first  input  segment  number 

source  node  for  the  first  in- 
put segment 

hot  output  segment  number 
=  -1  for  a  cold  stream 

destination  node  for  the  hot 

output  segment 

=  -1  for  a  cold  stream 

second  input  segment  number 

source  node  for  the  second 
input  segment 
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Card 


Variable      Card  Format    Description 


NARRAY(I,8) 

NARRAY(1,9) 

NARRAY(i,10) 
NARRAY(I,11) 


Optional 

Cards  (when   ICON (I) 

NSEGS  +  I  =  1,NSEG 

0) 

Optional  SEC0N(I,J) 

Cards  J  =  1,2 

(when  NSECS  I  =  1, NSEGS 
4   0) 


1114 


8F8.1 


cold  output  segment  number 
=  -1  tor  a  hot  stream 

destination  node  for  the  cold 

output  segment 

=  -1  for  a  hot  stream 

type  of  node 
=  5  for  a  mixer 

Stream  I.D.  for  the  cold  stream 

segments 

=  0  for  a  hot  stream 

Note:   the  order  for  the  input 
segments  can  be  arbitrary. 

This  vector  contains  the  segment 
numbers  for  which  both  the  lower 
and  upper  bounds  are  to  be 
specified  for  enthalpies. 

SEC0N(1,1)  contains  the  lower 
bounds  on  enthalpies  correspond- 
ing to  the  segment  specified  at 
the  same  1  location  in  ICON. 

SEC0N(I,2)  contains  upper  bounds. 

8  entries  per  card  with  all  the 
lower  bounds  specified  first. 
Start  with  a  new  card  for  upper 
bounds . 


Explanation  About  the  Output  From  EROS 
The  output  provides  the  following  information: 

1)  A  printout  of  the  input. 

2)  The  progress  of  optimization  (optional) . 

3)  Incidence  matrix.   The  rows  represent  functions  and  columns 
represent  variables. 

4)  Information  about  functions.   The  first  column  points  to  the  row 
of  the  incidence  matrix  representing  the  function.   The  second 
column  lists  the  node  associated  with  the  function.   The  third 
column  indicates  the  type  of  equation  involved. 
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TYPE  =  1  Heat  balance  for  an  exchanger,  heater  or 

cooler 

TYPE  =  2,3  Material  balances  for  a  splitter 

TYPE  =  4  Material  balance  for  a  mixer 

TYPE  =  5  Heat  balance  for  a  mixer 

TYPE  >  10,  TYPE  <  0   Type  is  an  inequality  constraint  code. 
(See  Table  1.) 

5)  Identification  of  unknown  variables.   The  first  column  stands  for 
the  segment  involved.   If  the  flow,  temperature,  or  the  split 
fraction  associated  with  this  segment  is  an  unknown  variable,  the 
entry  indicates  the  column  number  in  the  incidence  matrix.   Several 
variables  may  have  the  same  column  number  if  the  system  discovers 
they  must  be  equal  by  heat  or  material  balance.   If  the  variable 

is  known,  a  value  of  0  is  used. 

6)  Segment  Information  -  Final  values  at  optimal  point. 

7)  Precedence  order  list  -  Together  with  8,  indicates  how  the  equa- 
tions were  solved  during  the  last  optimization  step. 

8)  Decisions  and  tears.   If  a  variable  is  a  decision,  the  correspond- 
ing value  for  its  Function  and  MLIST  columns  are  -1  and  1  respec- 
tively.  If  a  tear  variable  is  involved,  the  function  it  is 
assigned  to  is  indicated  in  the  corresponding  column. 

MLIST  =  2  represents  explicit  assignment. 
MLIST  =  3  represents  explicit  assignment. 

9)  Objective  function  -  Final  value. 

10)  Exchanger  Areas  -  FinaL  values. 

11)  Information  about  time  and  number  of  iterations  required. 
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APPENDIX  B 
DESCRIPTION  OF  THE  PROGRAM 

An  abbreviated  block,  diagram  of  the  program  is  illustrated  in 
Figure  23.   Only  the  subroutines  performing  important  functions  are 
shown.   A  typical  run  of  EROS  takes  place  in  the  following  manner. 
The  main  program  makes  a  call  to  the  subroutine  EXEC.   EXEC  reads  in 
the  input  data.   Next,  the  routine  START  is  called  to  perform  an 
initialization  wherein  all  the  undefined  parameters  are  set  to  a  value 
of  zero.   Once  this  has  been  accomplished,  EXEC  performs  the  task  of 
identifying  all  the  unknown  variables  and  numbering  them  sequentially 
by  using  the  routines  FLOCTS,  HCNTS ,  KNOWN  and  CLEAR.   EXEC  issues 
directives  to  find  the  first  solution  procedure  for  the  system.   With 
the  aid  of  FMATRX,  an  incidence  matrix  is  created  with  rows  representing 
model  equations  and  columns  representing  the  unknown  variables.   PRECED 
uses  this  mati ix  to  find  a  solution  procedure.   At  this  point  the 
control  is  transferred  from  EXEC  to  USPONT,  a  subroutine  coordinating 
the  two  functions,  optimization  and  creation  of  all  subsequent  solution 
procedures . 

Every  time  a  new  solution  procedure  is  created,  USPONT  trans- 
fers control  to  the  optimization  routine  USOPT.   USOPT  makes  use  of 
the  complex  method  represented  by  OPT,  and  the  initialization  for  the 
optimization  is  performed  by  ICOPT.   OPT  hands  back  different  sets  of 
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Fig.  23.   Description  of  the  Program  EROS 
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decision  variable  values  to  USOPT.   For  each  set  of  decision  variable 
values,  model  equations  are  solved  and  the  objective  function  found. 
USOPT  invokes  the  use  of  ANALYZ  in  order  to  solve  the  model  equations. 
The  convergence  of  the  iteration  variables  in  ANALYZ  is  aided  by  the 
Harwell  Library  subroutine  VA05AD.   Constraint  violations  are  checked 
for  in  SUBANZ.   If  no  constraints  are  violated,  USOPT  requests  for 
the  total  system  cost  via  the  routine  COST.   If  constraint  violations 
occur,  control  is  transferred  back  to  USPONT.   A  decision  regarding 
which  inequality  constraints  are  to  be  present  as  equalities  in  the 
equation  set  is  made  in  POINT,  and  then  a  new  incidence  matrix  and 
consequently  a  new  solution  procedure  is  found  in  USPONT.   Optimiza- 
tion is  started  once  again,  using  USOPT. 

When  the  stopping  criteria  for  optimization  are  satisfied, 
control  is  transferred  to  EXEC  from  USOPT  through  USPONT.   An  output 
is  requested  by  a  call  to  PRINT  and  the  run  terminated. 

A  brief  description  of  all  the  subroutines  and  function 
subprograms  used  in  EROS  is  presented  next.   The  modules  are  presented 
in  an  alphabetical  order  and  unless  they  are  specified  explicitly  as 
function  subprograms,  they  represent  subroutines. 


ACYCLIC 

This  routine  helps  in  finding  a  solution  procedure.   It 
consists  of  the  acyclic  algorithm. 

ADDER 

When  a  new  solution  procedure  is  found,  constraints  in  the 
equation  set  corresponding  to  slack  variables  whose  values  are  zero, 
are  held  tight. 

AHEI P 

Using  the  area  of  a  heat-exchanger  zone,  and  the  user 
specified  constants,  cost  is  calculated  in  this  subroutine. 
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ANALYZ 

Groups  of  equations  that  can  be  solved  without  iterating  and 
calculations  involving  an  iteration  loop  are  identified  so  that  the 
model  equations  may  be  solved. 

Subroutines  used:   LOOP 
INIT1 
VA05AD 
SUBANZ 
FTLLX 

AREA 

Areas  are  calculated  for  different  phase  zones  and  subsequently 
the  cost  for  an  exchanger  is  evaluated. 

Subroutine  used:    AHELP 

BUBPT 

This  is  a  function  subprogram  that  discovers  the  bubble-point 
temperature  for  a  stream  segment. 

CALFUN 

This  routine  is  called  by  VA05AD  to  provide  the  value  of  the 
error  function. 

Subroutine  used:    SUBANZ 

CENTRE 

A  centroid  is  found  for  the  points  generated  by  the  complex 
method . 

CLEAR 

The  unknown  variables  are  numbered  sequentially. 

Subroutine  used:    CLEAR1 

CLEAR1 

An  aid  to  the  routine  CLEAR. 

C0LOC 

This  subroutine  is  used  to  find  the  value  of  a  slack  variable 
from  the  current  best  point. 

C0N0DE 

The  unknown  variable  involved  in  a  constraint  present  as  an 
equality  is  calculated  by  this  routine. 

Subroutines  used:   FILLUP 
ENTHAL 
ENSLAK 
TEMP 
COLOC 
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Function  Subprograms  Used:   DEWPT 

CONSTR 

A  check  is  made  to  see  if  any  constraints  are  violated.   If 
none  are  violated  the  slack  values  for  all  the  constraints  are  stored. 

Subroutines  used:   ENTHAL 
TEMP 
ENSLAK 

Function  Subprograms  used:   DEWPT 

BUBPT 

CONVRT 

In  this  routine,  an  array  that  provides  information  about 
slack  variables  associated  with  constraints  in  the  equation  set,  is 
created. 

COST 

Once  all  the  model  equations  are  solved,  the  costs  associated 
with  exchangers  and  the  utility  streams  are  added  together  in  this 
subroutine  to  provide  the  value  of  the  objective  function. 

Subroutine  used:    AREA 

DEWPT 

This  is  a  function  subprogram  that  discovers  the  dew-point 
temperature  for  a  stream  segment. 

ENBALH 

During  the  setting  up  of  an  incidence  matrix,  this  routine 
helps  in  creating  a  row  for  the  heat-balance  of  an  exchanger. 

ENSLAK 

This  routine  helps  in  treating  the  temperature  constraints 
in  terms  of  enthalpies 

Subroutines  used:   ENTHAL 
TEMP 

ENTHAL 

Used  for  calculating  the  enthalpy  and  vapor  fraction  of  a 
segment  given  its  temperature.   For  a  pure  component  enthalpy  is 
calculated  given  both  temperature  and  vapor  fraction. 

Function  Subprograms  used:   BUBPT 

DEWPT 

EQS 

Helps  relate  flows,  enthalpies  or  areas  for  reliability 
analysis. 

EXEC 

This  is  an  executive  routine  to  which  all  the  input  data  are 
provided.   Calls  are  made  to  several  subroutines  for  the  execution  of 
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the  program.   The  directive  for  providing  the  final  output  is  also 

issued  by  EXEC. 

Subroutines  used:   START 

FLOCTS 

HCNTS 

KNOWN 

CLEAR 

FMATRX 

SETUP 

PRECED 

LINK 

ANALYZ 

COST 

USPONT 

PRINT 

FILLUP 

Once  an  unknown  variable  is  calculated,  its  value  is  stored 

for  all  the  segments  it  occurs  in. 

FILLX 

The  value  of  the  decision  variable  is  stored  for  all  the 

segments  this  variable  occurs  in. 

Subroutine  used:    TEMP 

FLOCTS 

All  the  segments  having  a  common  flow-rate  are  identified 

and  characterized  by  a  unique  number  in  this  routine. 


Subroutine  used: 


EQS 


FMATRX 

An  incidence  matrix  is  created  in  which  the  rows  represent 
the  heat  and  material  balances  and  inequality  constraints  converted 
to  equalities.   The  columns  represent  the  unknown  variables. 

Subroutine  used:    SETUP 

FPASS 

This  routine  is  used  when  a  new  solution  procedure  is  to  be 
created  by  replacing  one  of  the  constraints  present  in  the  equation 
set  by  the  most  recently  violated  constraint. 


HCNTS 

All  the  segments  having  a  common  enthalpy  are  characterized 
by  a  unique  number  in  this  routine. 


Subroutine  used; 


EQS 
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HTEX 

Any  unknown  variable  in  the  heat-balance  equation  of  the 
heat-exchanger  is  calculated  by  this  subroutine.   The  error  function 
associated  with  the  heat-balance  equation  is  also  calculated. 

Subroutine  used:   FILLUP 
TEMP 

HTTC 

This  function  subprogram  is  used  to  calculate  the  overall 
heat-transfer  coefficient  based  on  the  film  heat-transfer  coefficients 
of  the  streams  involved. 

ICOPT 

Values  for  decision  variables  and  other  parameters  involved 
in  optimization  are  initialized. 

Subroutines  used:   C0L0C 
SORT 

INIT1 

Initializes  parameters  for  the  use  of  VA05AD. 

KNOWN 

In  this  routine  the  known  values  for  flows  and  enthalpies 
are  taken  into  account  to  help  discover  what  the  unknown  variables 
are. 

Subroutine  used:   ENTHAL 

LINK 

An  array  is  created  with  entries  that  keep  track  of  the 
segments  that  have  the  same  flow  or  enthalpy. 

LOOP 

Identifies  all  the  iteration  loops  for  ANALYZ . 

MB  HAD 

A  Harwell  Library  subroutine  called  by  VA05AD. 

MIXER 

The  unknown  variables  involved  in  the  heat  and  material 
balances  of  a  mixer  and  the  error  functions  are  calculated  by  this 
subroutine. 

Subroutines  used:   TEMP 

FILLUP 

OPT 

This  routine  consists  of  the  modified  complex  algorithm 
for  optimization. 

Subroutine  used:   CENTRE 
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POINT 

In  this  routine  a  decision  is  made  regarding  what  constraints 
are  to  be  incorporated  into  the  equation  set  to  obtain  a  new  solution 
procedure. 

PRECED 

A  solution  procedure  is  discovered  for  the  given  set  of 
equations  represented  in  the  matrix  form.   The  decisions  and  tear 
variables  are  established. 

Subroutine  used:   ACYCLIC 

PRINT 

Provides  an  output. 

RECOG 

An  array  is  created  in  which  information  concerning  the  un- 
known variables  (whether  they  represent  flows,  enthalpies,  or  split 
fractions)  is  stored. 

RELEAS 

Slack  variables  corresponding  to  constraints  in  the  equation 
set  are  identified  if  they  are  to  be  treated  as  search  coordinates. 

SAFE 

The  value  of  the  objective  function  corresponding  to  the 
best  current  point  is  stored. 

SETUP 

Creates  an  incidence  matrix  with  rows  represented  only  by 
the  heat  and  material  balance  equations. 

Subroutine  used:   ENBALH 

SORT 

A  routine  that  creates  an  array  containing  information  about 
slack  variables  that  are  to  be  used  as  search  coordinates. 

SPLITR 

Any  unknown  variables  in  the  material  balance  equations  of 
a  splitter,  and  the  error  functions  involved,  are  calculated  by  this 
subroutine. 

Subroutine  used:   FILLUP 

START 

Initializes  all  the  unspecified  parameters  at  the  start  of 
execution. 

Subroutine  used:   EQS 
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SUBANZ 

For  the  calculation  of  an  unknown  variable,  a  call  is  made 
to  the  unit  model  routines  associated  with  the  variable.   After  all 
the  calculations  have  been  completed,  CONSTR  is  called  to  check  for 
constraint  violations. 

Subroutines  used:   HTEX 

SPLITR 
MIXER 
C0N0DE 
CONSTR 

TEMP 

The  temperature  and  vapor  fraction  of  a  segment  are  calcu- 
lated given  its  enthalpy. 

USOPT 

Optimization  is  performed  using  the  strategy  followed  in 
this  study. 

Subroutines  called:   CONVRT 
ADDER 
RELEAS 
1C0PT 
OPT 

ANALYZ 
COST 

USPONT 

The  algorithm  for  finding  a  new  solution  procedure  is  co- 
ordinated with  the  optimization  strategy. 

Subroutines  used:   SETUP 
FMATRX 
PRECED 
LINK 
RECOG 
USOPT 
FPASS 
POINT 
PRINT 

VA05AD 

A  Harwell  Library  subroutine  that  performs  a  least-square 
fit  on  a  set  of  non-linear  equations. 

Subroutines  used:   MB11AD 

CALFUN 

In  addition,  a  system  subroutine  (CPUTIM)  on  IBM  360/67  ar  Carnegie- 
Mellon  University  was  used.   This  routine  provided  the  elapsed  CPU 
time  in  milliseconds.   It  was  called  by  the  following  subroutines: 

EXEC 

ANALYZ 

USPONT 
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